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Abstract: We describe an algorithm to construct an intrinsic Delaunay triangulation of a 
smooth closed submanifold of Euclidean space. Using results established in a companion paper 
on the stability of Delaunay triangulations on 5-generic point sets, we establish sampling criteria 
which ensure that the intrinsic Delaunay complex coincides with the restricted Delaunay complex 
and also with the recently introduced tangential Delaunay complex. The algorithm generates a 
point set that meets the required criteria while the tangential complex is being constructed. In 
this way the computation of geodesic distances is avoided, the runtime is only linearly dependent 
on the ambient dimension, and the Delaunay complexes are guaranteed to be triangulations of the 
manifold. 
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Resume : Nous montrons que, pour toute sous variete M suffisamment reguliere de l'espace 
euclidien et pour tout echantillon P de points de M qui satisfait une condition locale de delta- 
genericite et de epsilon-densite, P admet une triangulation de Delaunay intrinseque qui est egale 
a la triangulation de Delaunay restreinte a M et aussi au complexe de Delaunay tangent. Nous 
montrons egalement comment produire de tels ensembles de points. 
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1 Introduction 

This paper addresses the problem of constructing an intrinsic Delaunay triangulation of a smooth 
closed submanifold A4 C R w . We present an algorithm which generates a point set V C M. and 
a simplicial complex on V that is homeomorphic to A4 and has a connectivity determined by 
the Delaunay triangulation of V with respect to the intrinsic metric of M. 

For a submanifold of Euclidean space, the restricted Delaunay complex |ES97_, which is de- 
fined by the ambient metric restricted to the submanifold, was employed by Cheng et al. [CDR05 
as the basis for a triangulation. However, it was found that sampling density alone was insufficient 
to ensure a triangulation, and manipulations of the complex were employed. 

In an earlier work, Leibon and Letschcr [LL00_ announced sampling density conditions which 
would ensure that the Delaunay complex defined by the intrinsic metric of the manifold was a 
triangulation. In fact, as shown in Section |2.4.3| and Appendix [A] the stated result is incor- 
rect: sampling density alone is insufficient to guarantee an intrinsic Delaunay triangulation (see 



Theorem A. 3 1. Topological defects can arise when the vertices lie too close to a degenerate or 
"quasi-cospherical" configuration. 

Our interest in the intrinsic Delaunay complex stems from its close relationship with other 
Delaunay-like structures that have been proposed in the context of non-homogeneous metrics. For 
example, anisotropic Voronoi diagrams [LS03J and anisotropic Delaunay triangulations emerge 
as natural structures when we want to mesh a domain of K m while respecting a given metric 
tensor field. 

This paper builds over preliminary results on anisotropic Delaunay meshes [BWYll] and 
manifold reconstruction using the tangential Delaunay complex [BG11 . The central idea in both 
cases is to define Euclidean Delaunay triangulations locally and to glue these local triangulations 
together by removing inconsistencies between them. We view the inconsistencies as arising from 
instability in the Delaunay triangulations, and exploit the results of a companion paper [BDG^ 
to define sampling conditions under which these inconsistencies cannot arise. 

The algorithm is based on the tangential Delaunay complex |BG11| . and is an adaptation 
of a Delaunay refinement algorithm designed to avoid poorly shaped "sliver" simplices [Li03 
IBGlOj . The tangential Delaunay complex is defined with respect to local Delaunay triangulations 
restricted to the tangent spaces at sample points. We demonstrate that the algorithm produces 
sampling conditions such that the tangential Delaunay complex coincides with the restricted 
Delaunay complex and the intrinsic Delaunay complex. The refinement algorithm avoids the 
problem of slivers without the need to resort to a point weig hting strategy (CDE+00[ ICDR05 
IBG11| , which alters the definition of the restricted Delaunay complex. 

We present background and foundational material in Section[2] Then, in Section|3] we exploit 
results established in [BDG12 to demonstrate sampling conditions under which the intrinsic De- 
launay complex, the restricted Delaunay complex, and the tangential Delaunay complex coincide 
and are manifold. The algorithm itself is presented in Section|4] and the analysis of the algorithm 
is presented in Section [5] 



2 Background 

Within the context of the standard m-dimensional Euclidean space K m , when distances are 
determined by the standard norm, ||-||, we use the following conventions. The distance between a 
point p and a set X c M m , is the infimum of the distances between p and the points of X, and is 
denoted dRm(p, X). We refer to the distance between two points a and b as ||&— a\\ or dRm(a, 6) 
as convenient. A ball -Br™(c, r) = {x \ \\x — c\\ < r} is open, and Bgm(c, r) is its topological 
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closure. Generally, we denote the topological closure of a set X by X, the interior by int(X), 
and the boundary by dX. The convex hull is denoted conv(X), and the affine hull is aff(X). 

We will make use of other metrics besides the Euclidean one. A generic metric is denoted d, 
and the associated open and closed balls are B(c, r), and B(c, r). If a specific metric is intended, 
it will be indicated by a subscript, for example in Section [3] we introduce d_M, the intrinsic metric 
on a manifold A4, which has associated balls £?jh(c, r). 

If A is a k x j matrix, we denote its i th singular value by Si(A). We use the operator norm 
\\A\\ = Sl (A) =sup w=1 \\Ax\\. 

If U and V are vector subspaces of E m , with dim U < dim V, the angle between them is 
defined by 

sinZ([/, V) — sup \\u — 7ryu|| , 

u£U 

where is the orthogonal projection onto V. This is the largest principal angle between U and 
V. The angle between affine subspaces K and H is defined as the angle between the corresponding 
parallel vector subspaces. 

2.1 Sampling parameters and perturbations 

The structures of interest will be built from a finite set P C K m , which we consider to be a set of 
sample points. If D C R m is a bounded set, then P is an e-sample set for D if d-Rm(x, P) < e for 
all x £ D. We say that e is a sampling radius for D satisfied by P. If no domain D is specified, 
we say P is an e-sample set if d(x, P U dconv(P)) < e for all x £ conv(P). Equivalently, P is an 
e-sample set if it satisfies a sampling radius e for 

D e (P) = {x £ conv(P) | d Rm (a;,5conv(P)) > e}. 

In particular, if V is an e-sample set for U, and P = U n V, and conv(P) C U, then P is an 
e-sample set. 

A set P is X-sparse if or™ (p, q) > A for all p,q £ P. We usually assume that the sparsity of a 
e-sample set is proportional to e, thus: A = fioe. 

We consider a perturbation of the points P C M. m given by a function £ : P — ► ]R m . If ( 
is such that d^m.(p, £(p)) < A we say that £ is a p -perturbation. As a notational convenience, 
we frequently define P = C(P)j an d let p represent G P. We will only be considering 

p-perturbations where p is less than half the sparsity of P, so ( : P — > P is a bijection. 

Points in P which are not on the boundary of conv(P) are interior points of P. 

2.2 Simplices 

Given a set of j + 1 points {po, ■ • ■ >Pj} C P C K m , a (geometric) j-simplex a = [po, . . . ,pj] is 
defined by the convex hull: a — conv({po, • ■ ■ iPj})- The points pi are the vertices of a. Any 
subset {pi , . . • ,Pi k } of {po, . . . ,pj} defines a fc-simplex r which we call a /ace of er. We write 
t < (7 if r is a face of cr, and r < a if r is a proper face of c, i.e., if the vertices of r are a proper 
subset of the vertices of a. 

The boundary of cr, is the union of its proper faces: do- = [J T<a t. In general this is distinct 
from the topological boundary defined above, but we denote it with the same symbol. The 
interior of a is int(er) = a \ da. Again this is generally different from the topological interior. 
Other geometric properties of a include its diameter (the length of its longest edge), A (it), and 
the length of its shortest edge, L(a). If a is a 0-simplex, we define L(a) = A(cr) = 0. 
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For any vertex p € a, the face oppposite p is the face determined by the other vertices of a, 
and is denoted a p . If r is a j-simplex, and p is not a vertex of r, we may construct a (J + 1)- 
simplcx it = p * r, called the join of p and r. It is the simplex defined by p and the vertices of 
r, i.e., r = a p . 

Our definition of a simplex has made an important departure from standard convention: 
we do not demand that the vertices of a simplex be affinely independent. A ^'-simplex a is a 
degenerate simplex if dimaff(er) < j. If we wish to emphasise that a simplex is a j-simplex, we 
write j as a superscript: o~ J ; but this always refers to the combinatorial dimension of the simplex. 

If a is non-degenerate, then it has a circumcentre , C{a), which is the centre of the smallest 
circumscribing ball for a. The radius of this ball is the circumradius of cr, denoted R(a). A 
degenerate simplex may or may not have a circumcentre and circumradius. We write R(o~) to 
indicate that a simplex has a circumcentre. We will make use of the afiine space iV(cr) composed 
of the centres of the balls that circumscribe cr. We sometimes refer to a point c 6 N(o-) as a 
centre for a. The space N(cr) is orthogonal to aff (cr) and intersects it at the circumcentre of a. 
Its dimension is m — dim aff (cr). 

The altitude of p in a is D(p,a) = c?Rm (p, aff (cr p )). A poorly-shaped simplex can be char- 
acterized by the existence of a relatively small altitude. The thickness of a j-simplex a is the 
dimensionless quantity 



We say that cr is To-thick, if T(cr) > To. If a is To-thick, then so are all of its faces. Indeed if 
t < <7, then the smallest altitude in r cannot be smaller than that of cr, and also A(r) < A(cr). 

Although he worked with volumes rather than altitudes, Whitney |Whi57) p. 127] proved 
that the affinc hull of a thick simplex makes a small angle with any hyperplane which lies near 
all the vertices of the simplex. We can state this [BDG12 I Lemma 2.5] as: 

Lemma 2.1 (Whitney angle bound) Suppose cr is a j-simplex whose vertices all lie within 
a distance n from a fc-dimensional afiine space, H C M m , with k > j. Then 



2.2.1 Simplex perturbation 

We will make use of two results displaying the robustness of thick simplices with respect to small 
perturbations of their vertices. The first observation bounds the change in thickness itself under 
small perturbations: 

Lemma 2.2 (Thickness under perturbation) Let a = [po, . . . ,pj] and a = [po, . . . ,pj] be 
j-simplices such that — Pi\\ < p for all i € {0, . . . , j}. For any positive r) < 1, if 




if j = 
otherwise. 



sinZ(aff(cr),iJ) < 



2v 



T(a)A(aY 



(l- V )T(a) 2 L(a) 



(1) 



14 



then 



D(pi,a) > T]D(pi,a) 



for all i £ {0, . . . , j}. It follows that 



T(cr)A(cr) > 77T(cr)A(cr), 
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and 



T(or) > 1 - 



_2p_ 



r?T(a) > - V T(a). 



Proof Let p,q £ a with p, q the corresponding vertices of a. Let v — p — q and v = 
Define # = Z(v,aS(a p )) and = Z(v, aff (tip)). Since T(ct) < T(<7 p ), Whitney's Lemma 
us bound Z(aff(<7 p ), aff (<5p)) by the angle a defined by 

2 P 



2.1 



lets 



T(a)A(a) 



Also, by an elementary geometric argument, 



sin 7 



2 f > 



defines 7 as an upper bound on the angle between the lines generated by v and v. 
Thus we have 

D{p, a) = \\v\\sm6> i\\v\\ - 2p) sin(0 - a - 7) . 

Using the addition formula for sine together with the facts that for x, y <G [0, (1 — x) < cos 2; 
2sinx > x; and sinx + siny > sin(x + y), we get 



1 - 2 



Dip, d) > i\\v\\-2p) 
For convenience, define p 

Dip, a) > |H|(l-// i 

> (i-m) 

> (i-m) 

>(l--^- 

> jsT£>(p, cr) 



2p 



2p\\ Dip, a) 



T(a)A(a) || V || 



2 P > _?p_ > 2 P Then 
L(<t) - ||t)|| - A(»- lnCn 



1-21 



2p 



T(a)A(a) + || V || 



T(a) 



T(a) 



Dip, a) 
Dip, a) 



T(a) 
2p||u 



V 2p 



T(a) T(a) 2 



Dip, a) 



when p < 



T(^) 2 A(a) 
D(p,<7) 

(l-K)T(a) 2 



£>(p,<7) 



The condition on p is satisfied when p satisfies Inequality (JT|). 

The bound on Y(er)A(<7) follows immediately from the bounds on the Dip, a), and the bound 
on T(ct) itself follows from the observation that 



A(a) > A(cr) 



A(or) " A(tr) + 2p 
when p satisfies Inequality ([!]). 



> 1- 



2p 
A(a) 



> 1- 



T(a) 



2^ 6 
^7< 



□ 



We will also make use of a bound relating circumscribing balls of a simplex that undergoes a 
perturbation: 
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Lemma 2.3 (Circumscribing balls under perturbation) Let a = [po, . . . ,pj] and a = [po, ■ 
be j-simplices such that \\pi — pi\\ < p for all i £ {0, . . . , j}. Suppose B = _Brto(c, r), with r < e, 
is a circumscribing ball for a. If 

T(<7) 2 L(<7) 

P - 28 

then there is a circumscribing ball i? = B R d (c, f) for a with 

l |e - c|l< (w))" < 2 > 

and 

If, in addition, we have that po = po, then |f — r| < ||c — c||, and ([2]) serves also as a bound on 
\f — r\. 

Proof By the perturbation bounds, the distances between c and the vertices of a differ by no 
more than 2p. Also, ||c — pi\\ < e = e + p, and so by [BDG121 Lemma 4.3] we have 

d um (c,N(a))< -'<' 



T(or)A(or)- 



The bound on p allows us to apply Lemma 2.2 with K = |, so T(cr)A(cr) > ^T(cr)A(cr), and we 
obtain the bound on ||c — c|| with the observation that e < 2e. Indeed, p < e because L(a) < 2e. 

By the triangle inequality \f — r\ < ||po — Poll + 1 1 2 — c|| , and the stated bound on \r — r\ 
follows from the observation that x(er)A(o-) ^ 1 if J > 1- Under the assumption that p — Poi the 
bound on lie — ell also serves as a bound on If — r\. □ 



2.2.2 Flakes 

For algorithmic reasons, it is convenient to have a more structured constraint on simplex geometry 
than that provided by a simple thickness bound. A simplex that is not thick has a relatively 
small altitude, but we wish to exploit a family of bad simplices for which all the altitudes are 
relatively small. As shown by Lemma |2.7| below, the r -flakes form such a family. The flake 
parameter r is a positive real number smaller than one. 

Definition 2.4 (To-good simplices and To-flakes) A simplex a is To-good if T(a^) > T Q for 
all j-simplices a 3 < a. A simplex is To-bad if it is not To-good. A Tq- flake is a To-bad simplex 
in which all the proper faces are T -good. 

Observe that a flake must have dimension at least 2, since T(a 3 ) — 1 for j < 2. A flake that 
has an upper bound on the ratio of its circumradius to its shortest edge is called a sliver. The 
flakes we will be considering have no upper bound on their circumradius, and in fact they may 
be degenerate and not even have a circumradius. 

Ensuring that all simplices are To-good is the same as ensuring that there are no flakes. 
Indeed, if a is To-bad, then it has a j-face a 3 < a that is not T -thick. By considering such a 
face with minimal dimension we arrive at the following important observation: 

Lemma 2.5 A simplex is To-bad if and only if it has a face that is a To-flake. 

We obtain an upper bound on the altitudes of a To-flake through a consideration of dihedral 
angles. In particular, we observe the following general relationship between simplex altitudes: 
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Lemma 2.6 If a is a j-simplex with j > 2, then for any two vertices p,q € a, the dihedral angle 
between o~ p and <j q defines an equality between ratios of altitudes: 

smZ aff cr p ,aff (cr q )) = — = — -. 

D{p,a q ) D(q,(Xp) 

Proof Let a pq = a p C\a q , and let p* be the projection of p into aS(a pq ). Taking p* as the origin, 
we see that ^ p v * has the maximal distance to aff(cTp) out of all the unit vectors in aff(cr g ), and 

this distance is gfe^j ■ By definition this is the sine of the angle between aff(<7 p ) and aS(a q ). 
A symmetric argument is carried out with q to obtain the result. □ 

We arrive at the following important observation about flake simplices: 

Lemma 2.7 (Flakes have small altitudes) If a fc-simplex a is a r -flake, then for every ver- 
tex p € a, the altitude satisfies the bound 

fcA(q) 2 r 

D{p > a)< (k-i)L(*y 

Proof Recalling Lemma [2. 6| we have 

DM = D{q $ D(p : a *\ 

D{q,a p ) 

and taking q to be a vertex with minimal altitude, we have 

D(q,a) = kT(a)A(a) < kT%A(a), 

and 

D(q, a p ) > (k - l)T(a p )A(a p ) > (k - l)^ 1 L(o-) , 

and 

D(p,a g )< A(a q )< A{a), 
and the bound is obtained. □ 

2.3 Complexes 

Given a finite set P, an abstract simplicial complex is a set of subsets K C 2 P such that if 
(7 € A", then every subset of a is also in A. The Delaunay complexes we study are abstract 
simplicial complexes, but their simplices carry a canonical geometry induced from the inclusion 
map 6 : P ^ E TO . (We assume l is injective on P, and so do not distinguish between P and 
6(P).) To each abstract simplex a £ A, we have an associated geometric simplex conv(6(<r)), and 
normally when we write a G A, we are referring to this geometric object. Occasionally, when it 
is convenient to emphasise a distinction, we will write b(a) instead of a. 

Thus we view such a A as a set of simplices in W n , and we refer to it as a complex, but it is not 
generally a (geometric) simplicial complex. A geometric simplicial complex is a finite collection 
G of non-degenerate simplices in M. N such that if a £ G, then all of the faces of a also belong 
to G, and if a, a 6 G and r = er n <7^0, then r < ct and t < a. An abstract simplicial complex 
is defined from a geometric simplicial complex in an obvious way. A geometric realization of 
an abstract simplicial complex A is a geometric simplicial complex whose associated abstract 
simplicial complex may be identified with A. A geometric realization always exists for any 
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complex. Details can be found in algebraic topology textbooks; the book by Munkres (Mun84 
for example. 

The carrier of an abstract complex A is the underlying topological space |A|, associated with 
a geometric realization of A . Thus if G is a geometric realization of A, then \K\ = {J aeG o. For 
our complexes, the inclusion map i induces a continous map i : \K\ — > R m , defined by barycentric 
interpolation on each simplex. If this map is injective, we say that A is embedded. In this case 
l also defines a geometric realization of A, and we may identify the carrier of A with the image 
of L. 

A subset A' C A is a subcomplex of A if it is also a complex. The star of a subcomplex 
K' C A is the subcomplex generated by the simplices incident to A' . I.e., it is all the simpliccs 
that share a face with a simplex of if', plus all the faces of such simpliccs. This is a departure 
from a common usage of this same term in the topology literature. The star of A' is denoted 
star(A') when there is no risk of ambiguity, otherwise we also specify the parent complex, as in 
star 

A triangulation of P C M. m is an embedded complex A with vertices P such that \K\ = 
conv(P). 

Definition 2.8 (Triangulation at a point) A complex K is a, triangulation at p £ R m if: 

• p is a vertex of A. 

• star(p) is embedded. 

• p lies in int(|star(p)|). 

• For all t £ A, and a £ star(p), if int(r) n a ^ 0, then r G star(p). 

A complex if is a j -manifold complex if the star of every vertex is isomorphic to the star of a 
triangulation of W . 

If a is a simplex with vertices in P, then any map £ : P — > P C M m defines a simplex C( fJ ) 
whose vertices in P are the images of vertices of it. If K is a complex on P, and K is a complex 
on P, then Q induces a simplicial map K — > K if C( cr ) £ ^ f° r every a £ K. We denote this map 
by the same symbol, £. We are interested in the case when £ is an isomorphism, which means it 
establishes a bijection between K and K . We then say that K and if are isomorphic, and write 
C - . P 

if = A , or K = K if we wish to emphasise that the correspondence is given by £. 
We use the following local version of a standard result [BDG12. Lemma 2.7]: 

Lemma 2.9 Suppose A is a complex with vertices P C K m , and K a complex with vertices 
P C K™. Suppose also that A is a triangulation at p £ P, and that ( : P — » P induces an injective 
simplicial map star(p) — >■ star(£(p)). If A is a triangulation at £(p), then 

C(star(p)) = star(C(p)). 

2.4 The Delaunay complex 

An empty ball is one that contains no point from P. 

Definition 2.10 (Delaunay complex) A Delaunay ball is a maximal empty ball. Specifically, 
B = i?K m {x, r) is a Delaunay ball if any empty ball centred at x is contained in B. A simplex a 
is a Delaunay simplex, if there exists some Delaunay ball B such that the vertices of a belong 
to dB n P. The Delaunay complex is the set of Delaunay simplices, and is denoted Dcl(P). 

The Delaunay complex has the combinatorial structure of an abstract simplicial complex, but 
Del(P) is embedded only when P satisfies appropriate genericity requirements |BDG12| . 
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2.4.1 Protection 

A Delaunay simplex a is 8-protected if it has a Delaunay ball B such that d^ [q, BE) > 8 for all 
q E P\cr. We say that B is a (^-protected Delaunay ball for a. If t < a. then B is also a Delaunay 
ball for r, but it cannot be a (5-protected Delaunay ball for r. We say that a is protected to mean 
that it is (5-protected for some unspecified <5 > 0. 

Definition 2.11 (^-generic) A finite set of points P C K m is S-generic if all the Delaunay 
m-simplices are (5-protected. The set P is simply generic if it is <5-generic for some unspecified 
6 > 0. 

We have previously demonstrated [BDG12| that (5-generic point sets impart a quantifiable 
stability on the Delaunay complex. In Section [3] we review the main stability result and develop 
it to define the sampling conditions that will be met by the algorithm that we introduce in 
Section |4] 

2.4.2 The Delaunay complex in other metrics 

We will also consider the Delaunay complex defined with respect to a metric d on M. m which 
differs from the Euclidean one. Specifically if P C U C M m and d : U x U — > E is a metric, then 
we define the Delaunay complex Delrf(P) with respect to the metric d. 

The definitions are exactly analogous to the Euclidean case: A Delaunay ball is a maximal 
empty ball B(x,r) in the metric d. The resulting Delaunay complex Del^P) consists of all the 
simplices which are circumscribed by a Delaunay ball with respect to the metric d. The simplices 
of Delrf(P) are, possibly degenerate, geometric simplices in M. m . As for Del(P), Del^P) has the 
combinatorial structure of an abstract simplicial complex, but unlike Del(P), Del^P) may fail 
to be embedded even when there are no degenerate simplices. 

2.4.3 Obtaining Delaunay triangulations in other metrics 

Delaunay |Del34 showed that if P C M. m is generic, then Del(P) is a triangulation. Point 
sets that are not generic are often dismissed in theoretical work, because an arbitrarily small 
perturbation of the points can be made which will yield a generic point set. Thus in the sense 
of the standard measure in the configuration space R mx l p l, almost all point sets will yield a 
Delaunay triangulation. However, when the metric is no longer Euclidean, this is no longer true. 

In contrast to the purely Euclidean case, topological problems arise in point sets that are 
"near degenerate" , i.e., point sets that are not (5-generic for a sufficiently large 5. How large <5 
needs to be depends on how much the metric differs from the Euclidean one. Indeed, this was 
the initial motivation for the introduction of (5-generic point sets |BDG12| , which are central to 
the results presented in this paper. 

As we show here with a qualitative argument, the problem can be viewed as arising from 
the fact that when m is greater than two, the intersection of two metric spheres is not uniquely 
specified by m points. We demonstrate the issue in the context of Delaunay balls. The problem 
is developed quantitatively in terms of the Voronoi diagram in Appendix [A] 

We work exclusively on a three dimensional domain, and we are not concerned with "boundary 
conditions"; we are looking at a coordinate patch on a densely sampled compact 3-manifold. 

One core ingredient in Delaunay's triangulation result |Del34| is that any triangle r is the 
face of exactly two tetrahedra. This follows from the observation that a triangle has a unique 
circumcircle, and that any circumscribing sphere for r must include this circle. The affine hull 
of t cuts space into two components, and if r € Del(P), then it will have an empty circumsphere 
centred at a point c on the line through the circumcentre and orthogonal to aff(r). The point c 
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Figure 1: In three dimensions, three closed geodesic balls can all touch three points, u,v,p, on 
their boundary and yet no one of them is contained in the union of the other two. 

is contained on an interval on this line which contains all the empty spheres for t. The endpoints 
of the interval are the circumcentres of the two tetrahedra that share t as a face. 

The argument hinges on the assumption that the points are in general position, and the 
uniqueness of the circumcircle for r. If there were a fourth vertex lying on that circumcircle, 
then there would be three tetrahedra that have t as a face, but this configuration would violate 
the assumption of general position. 

Now if we allow the metric to deviate from the Euclidean one, no matter how slightly, the 
guarantee of a well defined unique circumcircle for r is lost. In particular, If three spheres Si, 
S2 and S3 all circumscribe t, their pairwise intersections will be different in general. I.e., 

Si n s 3 ^ s 2 n s 3 . 

Although these intersections may be topological circles that are "arbitrarily close" assuming the 
deviation of the metric from the Euclidean one is small enough, "arbitrarily close" is not good 
enough when the only genericity assumption allows configurations that are arbitrarily bad. 

An attempt to illustrate the problem is given in Figure [I] where r = {u,v,p}. Here, the 
sphere S3 would be contained inside the spheres Si and S2 if the metric were Euclidean, but 
any aberration in the metric may leave a part of S3 exposed to the outside. This means that in 
principle another sample point w could lie on S3, while Si and 5*2 remain empty. Thus there are 
three tetrahedra that share r as a face. 

The essential difference between dimension 2 and the higher dimensions can be observed by 
examining the topological intersection properties of spheres. Specifically, two (m — l)-spheres 
intersect transversely in an (m — 2)-sphere. For a non-Euclidean metric, even if this property 
holds for sufficently small geodesic spheres, only in dimension two is the sphere of intersection 
of the Delaunay spheres of two adjacent 771-simplices uniquely determined by the vertices of the 
shared (m — l)-simplex. See Figure IT] 
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2.4.4 The Voronoi diagram 

We will occasionally make reference to the Voronoi diagram, which is a structure dual to the 
Delaunay complex. It offers an alternative way to interpret observations made with respect to 
the Delaunay complex. 

The Voronoi cell associated with p € P with respect to the metric d : U x U — > K is given by 

V d (p) = {x eU\ d(x,p) < d(x,q) for all q E P}. 

More generally, a Voronoi face is the intersection of a set of Voronoi cells: given {po, . . . ,pk} C P, 
let a denote the corresponding abstract simplex. We define the associated Voronoi face as 

k 

Vd{<r) = f| v ^)- 

It follows that a is a Delaunay simplex if and only if Vd(cr) ^ 0. In this case, every point in 
Vd(cr) is the centre of a Delaunay ball for a. Thus every Voronoi face corresponds to a Delaunay 
simplex. The Voronoi cells give a decomposition of U, denoted Vor^P), called the Voronoi 
diagram. Our definition of the Delaunay complex of P corresponds to the nerve of the Voronoi 
diagram. 



3 Equating Delaunay structures 

We now turn to the task of triangulating Ai, a smooth, compact m-manifold, without boundaries 



embedded in M. N . In this section we demonstrate our main structural result, Theorem 3.5 which 
is stated at the end of Section |3.1| It says that the complex constructed by the algorithm we 
describe in Section [4] is in fact an intrinsic Delaunay triangulation of the manifold, which we 
introduce next. 



3.1 Delaunay structures on manifolds 

The restricted Delaunay complex is the Delaunay complex Del^N^ (V) obtained when distances 
on the manifold are measured with the metric d-^N\ M . This is the Euclidean metric of the ambient 
space, restricted to the submanifold Ai. In other words, d^N^ M (x,y) — d^N(x,y). We use this 
notation to avoid ambiguities in conjunction with the local Euclidean metrics discussed below. 
The Delaunay complex DelgN^ (V) is a substructure of Del^^T 3 ). 

Alternatively, distances on the manifold may be measured with the intrinsic metric of 
the manifold. This metric defines the distance between x and y as the infimum of the lengths of 
the paths on Ai which connect x and y. Since the length of a path on Ai is defined as its length 
as a curve in WL N , this metric is also induced from cZrat. The intrinsic Delaunay complex is the 
Delaunay structure DeDvf('P) associated with this metric. 

Although neither of these metrics are Euclidean, the idea is that locally, in a small neigh- 
bourhood of any point, these metrics may be well approximated by d^m. Then, if the sampling 
satisfies appropriate (5-generic and e-dense criteria in these local Euclidean metrics, the global 
Delaunay complex in the metric of the manifold will coincide locally with a Euclidean Delaunay 
triangulation, and we can thus guarantee a manifold complex. 



3.1.1 Local Euclidean metrics 

A local coordinate chart at a point p G Ai, is a pair [W, 4> p ), where W C M. is an open 
neighbourhood of p, and (f> p : W — > U = <j) p (W) C K m is a homeomorphism onto its image, with 
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4>p{p) = 0. A local coordinate chart allows us to pull back the Euclidean metric to W. For all 
x, y £ W, the metric d,p p (x,y) = d^™ (4> p (x), (j> p (y)) is a a local Euclidean metric for p on W. 
This metric depends upon the choice of <f> p ; there are different ways to impose a Euclidean metric 
on W. 

It is convenient to take the reciprocal point of view, and work with a local parameterization 
at a point p £ M- This is a pair (U, ip p ), such that U C M. m , and (W, is a local coordinate 
chart for p, where W = ip p (U). We can then use %p p to pull back the metric of the manifold to 
U, and to simplify the notation we write dj^(x,y) for x,y £ U, where it is to be understood 
that this means dM('4 , p(x),^j p (y)), and likewise for d-^N^ M (x,y). Indeed, once W and U have 
been coupled together by a homeomorphism, we can transfer the metrics between them and the 
distinction becomes only one of perspective; the standard metric dgm on U is a local Euclidean 
metric for p. 

We wish to generate a sample set P C M that will allow us to exploit the stability results 
for Dclaunay triangulations [BDG12 . We consider the stability of a Delaunay triangulation in a 
local Euclidean metric. The following definition is convenient when stating the stability results: 

Definition 3.1 (Secure simplex) A simplex a £ Del(P) is secure if it is a (5-protected m- 
simplex that is T -thick and satisfies R{u) < e and L(o~) > n^e. 

We will make reference to the following result [BDG12, Theorem 4.17]: 

Theorem 3.2 (Metric stability assuming thickness) Suppose conv(P) CUc M m and the 

metric d : U x U — > K is such that \d(x,y) — du™(x,y)\ < p for all x,y £ U. Suppose also that 
Pj C P is such that every m-simplex a £ star(Pj; Del(P)) is secure and satisfies dgm (p, dU) > 2e 
for every vertex p £ a . If 

H ~ 36 ' 

then 

star(Pj; Del d (P)) = star(Pj; Del(P)). 



In our context the point set P used in Theorem |3.2| will come from a larger point set V, such 
that P = W n V . We will write P\y in order to emphasise this dependence on W. We want to 
ensure that 

star(P j; DeLj(P^)) = star(P, 7 ; Del d (T)). (3) 

This requirement is attained by demanding that V satisfy a sampling radius of e with respect 
to the metric c?». Since c£gm(a;, y) < d_M{x,y) for all x,y £ U = W, by our particular choice of 
ip p , we will have that P\y is an e-sample set with respect to the metric dm.™. We ensure that U 
is large enough so that d]gt™(p,dU) > 4e for all p £ Pj. It then follows that R(a) < e for any 
simplex a £ star(Pj; Del(Pw)), because Pw is an e-sample set [BDG12, Lemma 3.6], and thus 
d]Rm(g, 3C/) > 2e for any q £ a. It follows that d_\4(l,dU) > 2e as well, and thus the sampling 
radius on V ensures that Equation ^ is satisfied. For our purposes Pj will consist of a single 
point p, and the sampling radius e is constrained by the requirement that U be small enough 



that the metric distortion introduced by ip p meets the requirements of Theorem 3.2 



3.1.2 The tangential Delaunay complex 

The algorithm we describe in Section [4] is a variation of the algorithm described by Boissonnat 
and Ghosh [BG10 . This algorithm builds the tangential Delaunay complex, which we denote by 
Delrjvi('P). This is not a Delaunay complex as we have defined them, since it cannot be defined 
by the Delaunay empty ball criteria with respect to any single metric. However, it is a Delaunay- 
type structure, and as with DelgNi^ (V), the tangential Delaunay complex is a substructure of 
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DeljjN (V) . We will demonstrate sampling conditions which ensure that DelrMi'P) = Del^('P) = 
DeWi^P). 

Definition 3.3 (Tangential Delaunay complex) The tangential Delaunay complex for V C 
Ai C is defined by the criterion that a g Deljvvi^) ^ ^ has an empty circumscribing ball 
B k n (c, r) such that c G T p A4 for some vertex pGu. 

We define some local complexes to facilitate discussions of the tangential Delaunay complex. 
For all p e V, let 

K(p) = {a\V R N(a)nT p M^<6}, 

and define 

star(p) = star(p; K (p)). (4) 

Then the tangential Delaunay complex is the union of the complexes star(p) for all p G V . 

Boissonnat et al. |BG11[ Lemma 2.3] showed that V^n (V)PiT p Ai is equal to the m-dimensional 
weighted Voronoi diagram of V C T p M., where V is the orthogonal projection of V onto T p M. 
and the squared weight of a point p\ 6 V' is —\\pi — Pi\\ 2 - Therefore, K(p) is isomorphic to a 
dual complex (the nerve) of the fc-dimensional weighted Voronoi diagram of V ■ 

3.1.3 Power protection 

The algorithm introduced in Section |4.2| will ensure that for every simplex a in the tangential 
Delaunay complex, and every vertex p G cr, there is a Delaunay ball for a that is centred on 
T p A4 and is protected in the following sense: 

Definition 3.4 (Power protection) A simplex a with Delaunay ball B r n(C,R) is 5 2 -power- 
protected if d R N (C, q) 2 - R 2 > 5 2 for all q G V \ a. 

Observe that, if C $ M, the ball i? K « (C, R) is not an object that can be described by the metric 
d-R«\ M - In the context of the tangential Delaunay complex we use power-protection rather than 
the protection described in Section |2.4.1| because working with squared distances is convenient 
when we consider the Delaunay complex restricted to an affine subspace. 

3.1.4 Main structural result 

The rest of Section [3] is devoted to the proof of Theorem |3.5| below. It says that for the point set 
generated by our algorithm, the tangential Delaunay complex is isomorphic with the intrinsic 
Delaunay complex of M.. It then follows, from a previously published result [BG10, Theorem 5.1], 
that the intrinsic Delaunay complex is in fact homeormorphic to A4; it is an intrinsic Delaunay 
triangulation. 

Thus we obtain a partial recovery of the kind of results attempted by Leibon and Letscher [LLOO . 
Our sampling conditions, and our algorithm (existence proof) rely on the embedding of M. in 
R^; we leave purely intrinsic sampling conditions for future work. 

Theorem 3.5 (Intrinsic Delaunay triangulation) Suppose V C A4 is (jioe)-spa,ise with re- 
spect to g?kn, an d every m-simplex a G DelxMi'P) is To-thick, and has, for every vertex p G <r, a 
5 2 -power-protected empty ball of radius less than e centred on T p A4, with 5 > Sop-o 6 - If ^oAo — \i 
and 

f^^ 2 rch(X) 
€ ~ 1.5 x 10 6 ' 

then 

Vel TM {V) = De\ M N lM (V) = Bel M (V), 



Inria 



Constructing intrinsic DTs 



17 



and for e sufficiently small, these will be homeomorphic to A4: 

\De\ M (V)\=M. 

3.2 Choice of local Euclidean metric 

A local parameterization at p £ M. will be constructed with the aid of the orthogonal projection 

ir p :R N ^T p M, (5) 
restricted to M.. As shown in Lemma B.4 Niyogi et al. [NSW08, Lemma 5.4] demonstrated that 
if r < rch ^ ' , then ir p is a diffeomorphism from W — B^Ni M (p,r) onto its image U C T p M.. We 
will identify T p M. with M™, and define the homeomorphism 

1> P = n P \w :U—>W. (6) 

Using ip p to pull back the metrics d_M and d^N\ M to M. m , we can view them as perturbations of 
g?k™ • The magnitude of the perturbation is governed by the radius of the ball used to define W. 



Definition 3.6 We call a neighbourhood W of p £ M. admissible if W C B R N\ M (p,r), with 

< rch(Al) 
' — 100 • 

In all that follows, any mention of a local Euclidean metric refers to the one defined by tt p 
restricted to an admissible neighbourhood. The requirement r < ^j?^ is simply a convenient 
bound that yields a small integer constant in the perturbation bound of the following lemma, 
and does not constrain subsequent results. The bound could be relaxed to r < Ich ^- A/| ) at the 
expense of a weaker bound on the perturbation. 

Lemma 3.7 (Metric distortion) Suppose (U,ip p ) is a local parameterisation at p £ W C A4 
with W = ip p {U). UWC B m n 1m (p, r), with r < ££ ^ 11 , then for all x,y£U, 

23r 2 

\ M (x,y) -<hm (x, y) | < ^(x.y) -d R m(x,j/)| < ^.^ ■ 

Proof Let u, w € VU C -Brn^ (p, r), and let be the angle between the line segments [it, u] and 
[TT p (u),Tr p (v)], 0\ the angle between [u,v] and T U AL and 02 the angle between TpA4 and T U A4 
Thus < 0\ +02, and c?R>n (7r p (u), 7r p (w)) = d R « (it, w)cos 0. Defining 77 — ^f^y, Lemma 
yields 

d M (u,v) < d R N(u,v) (1 + 4n) , (7) 

and so 

<2m(u,v) cos6» 
d Rm (7Tp(u),7r p (w)) > . 



B.5 



Using Lemma B.l we find sin#i < 77, and Lemma [B.3| yields sin 02 < 677. Therefore, since 
sin6> < sin 0i + sin 6*2, we have cos = (1 — sin 2 6*) 1 / 2 > 1 — sin > 1 — 7rj and we get 



d M -m(n p (u),Tr p (v)) > d M (u,v) 



1 - 7r) 



1 + 

>d M (u,v)(l-7r])(l-4r]) 
> d M (u,v)(l- II77). 

Using Equation we find d M (u,v) < 2 U ^ : , so d R ™ (n p {u), tt p (v)) > d M (u,v) - 23 ^fbw) > and 
the result follows since g?jh(tj,?j) > d^N\ M (u,v) > d Rm (ir p (u), ir p (v)). □ 
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Our sampling radius is constrained by the size of a Euclidean ball that can be contained in 
an admissible neighbourhood. The following lemma gives a convenient expression for this: 

Lemma 3.8 If 1 < a < 10 4 and ae < rcl ^> , and U = B Mm (p, (a - l)e), then i/j p (U) =W C 



Proof Using Lemma B.l we have that B-gmfar) C w p (B^N\ M (p,ae)) if 

,2,2 x 2 / / x 2N 



2 ^ 2 2 

r < a e 



GrchCM)) ^f 1 (2rch(M)) 



< a 2 e 2 1 - 



200 / 



2N 



Thus we require r < J 2 °^ Q2 1 ae, which is satisfied by r = (a — l)e if a < 79999. □ 



Lemmas |3.7| and |3.8| lead to a sampling radius which allows us to employ Theorem |3.2| and 
so obtain an equivalence between Delaunay structures: 

Proposition 3.9 (Equating Delaunay complexes) Suppose V C M. is an e-sample set with 
respect to d^N\ M , and that for every p £ V , in the local Euclidean metric on W — B r n\ m (p, 5e), 
every m-simplex in star(p; Del(Pw)) is secure, where Pw = V D W, and S = v^e. If 

T fi iy ich(M) 
€ ~ 20700 

then 

star(p; Del(P w )) = star(p; Del R jv|^ (Pw)) = star(p; Del M (P w ))- (8) 

Thus 

Be\ miM (V) = Bcl M (V), 

and they arc manifold complexes. 



Proof As usual, let U = n p (W). Then by Lemma 3.8 -Br™(p, 4e) C U, and thus d^m.(q,dU) > 
2e for any vertex g of a simplex in star(p; Del(Pvy))- Thus Lemma 3.7 allows us to apply 
Theorem |3.2| provided 

23a 2 e 2 T Mo^o£ 



rch(TW) - 36 ' 

when a = 5, and we obtain the required bound on e. Thus the star of every vertex in Deljn('P) is 
equal to the star of that point in the local Euclidean metric, and likewise for ~Del^N\ M (V). The 
claim follows since a £ DeLvi (P) if an d only if it is in the local Euclidean Delaunay triangulation 
of every one of its vertices, and likewise for the simpliccs in DelgN^ (V). □ 



3.3 The protected tangential complex 

We obtain Theorem |3.5| by means of Theorem |3.2| via the observation that power protection of 
the ambient Delaunay balls translates into protection in the local Euclidean metrics. We must 
distinguish between the geometry of a simplex defined with respect to the Euclidean metric d R N 
of the ambient space, as opposed to a local Euclidean metric or™ . In general, we use a tilda to 
indicate simplices in the ambient space, and their properties. 
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Lemma 3.10 (Protection under projection) Suppose V C M. and that a € DelRjv('P) is an 
Tg-thick rn-simplex. with L(p) > /Etoe and B k n(C,R) is a 5 2 -power-protected empty ball for cx, 
with respect to the metric d^N , where <5 2 > ^/i^e 2 . Suppose also that C € T p Ai, for some vertex 

If i? < e, with 

f§A^ 2 rch(X) 
" 512 ' lj 

then a = TT p (a) has a (S-protected Delaunay ball Bgn.(c, r) with respect to the local Euclidean 
metric d R ™ for p on any admissible neighbourhood W that contains B-^n\ m (p, 3e), and (5 = ^ e. 
with 

-o= d f t) . (10) 

Proof We first find a bound for d$m(C, c) and r. Let a = \po, ■ ■ ■ , p m ], and er = [po, ■ • ■ ,Pm] so 
that 7T p (pi) = and p = po = Po- We will first show that, near C, there is a circumcentre c for 
tr in the metric tigm. For any 6 tr, d R w(p,pi) < 2i?, and so by Lemma B.l we have 

, N 2i? 2 2e 2 



In order to apply Lemma 



2.3 



rch(X) rch(M) 

2e 2 ^ ^nAo e 

we require -j-™ < -"g-, or 



e < 



rch(A^) — 28 

f 2 A rch(X) 



5G 



which is satisfied by Equatio n (|9| . Since aff(er) = T p .M, the circumcentre c € T p .A/f is the closest 
point in N(a) to C, Lemma 2.3 yields 

\R-r\<d Rm (C, c) = d m (C, c) < ^- 16< " 



T Q fi Q rch(M) 



Now we seek a lower bound on the protection of B-gm(c, r). Suppose q E V \ a. We wish to 
establish a lower bound on dR™(c, q) — r, where q — Tr p (q). We may assume that <iRjv(C, q) < 3e, 
since otherwise q will lie outside of our region of interest . 

z = 2rcSDwr) ^ e ^ ne u PP er bound on d^N(q,q) given by Lemma B.l Then dR m (C,q) 2 > 
d RN (C, q) 2 -z 2 > R 2 + S 2 - z 2 . Thus 



d m m(C,?)-i?> - > — , 

d Rm (C, g) + R 4e 

since i?<e. Then d Rm (c, q) - r > {d M m (C, q) - d R m(C, c)) - (R+ \R - r\) > ^f- -2d ffi m(C,c). 
Putting this together, using <5 2 > 5 2 //Qe 2 , we get 

j , , Ara -a 81e2 32e 

^ (C ' 9) " r > U V ° " I^MF - foAorch^) 

In order to simplify away the final term, we demand 

32e l^a 



T /lorch(X) " 16 
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which is satisfied by Equation Q. Under this constraint, the central term is also seen to be less 

1 



than Yg^Q/ig, and we obtain 



5 > ^ Ao e 



□ 



Proposition | 3.9| re quires a thickness T and shortest edge bound /j^e for the simplex a C 
M. m , but Lemma 3.10 is expressed in terms of the corresponding quantities To and ^.qc for the 
corresponding simplex a C M. N . 

Lemma 3.11 (Simplex distortion under projection) Let a e Del R w(7- > ) be an m-simplcx 
as described in Lemma 3.10 and let a = ffpft) be its projection in the local Euclidean metric 
for p on any admissible neighbourhood that contains B^n\ m (j>, 2e), where p is a vertex of a. If 
e satisfies Equation and ^oAo — \i then L(a) > floe, where 

20, 



and T(er) > To, where 



Ik) 



T 



21 



MO: 



49 T °' 



Proof Since ir p (B R N\ M (p,2e)) C i3 R ™(p, 2e), it is sufficient to apply the Metric distortion 
lemma [3~8l with a = 3. 

For the shortest edge length, we find 



L(a) > L{a) - 
= Mo U 

> Ao ( i 



20, 
> ^Moe. 



3 2 x 23e 2 
rch(M) 

512 



For the thickness bound, in order to apply Lemma 2.2 using rj = (1 — 77), we require 



207T ^ < ffCfa 



512 



14 



which is satisfied if we choose 
Then Lemma [2 . 2| yields 



T(a)>^(l-6^)T(a)>^T(a). 



□ 



We can now express the sampling conditions in terms of the output parameters of the tan- 
gential complex algorithm, and this allows us to apply Proposition 3.9 and obtain our main 
structural result: 
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Proof of Theorem \3.5\ We first translate the sampling requirements of Proposition |3.9| in terms 
of properties of simplices in the ambient metric d K w . Using Lemma |3.11[ together with Equa- 



tion (10 1, the upper bound on the sampling radius demanded by Proposition 3.9 becomes 



20 x 6fo<5g/Zgrch(M) 
~ 21 x 49 x 8 x 20700 ' 



We obtain the stated sampling radius bound after multiplying by T in order to ensure that the 
demand of Equation (|9| is also met. Thus the stated sampling radius satisfies the requirements 
of both Lemma |3.10| and Proposition |3.9| 

The fact that the structures are isomorphic follows from the fact that they are all lo- 
cally isomorphic to the Delaunay triangulation in the local Euclidean metric. To see that 
star(p; DelxM ('P)) — star(p; Del(Pv^)), observe that Lemma |3.10 implies that there is an in- 



jective simplicial map star(p; Delrx (V)) — > star(p; Del(Pw))- The isomorphism is established 
by Lemma |2.9| once it is established that star(p; DelyjuCP)) is a triangulation at p. In fact 
star(p; DelrM (T 7 )) is isomorphic to the star of p in a regular triangulation of the projected 
points Pw', it is a weighted Delaunay triangulation [BG10J, and with our choice of W, the point 
p is an interior point in this triangulation [BG10I Lemma 2.7(1)]. Thus star(p; DelrM i'P)) is a 
triangulation at p, and it follows that 

star(p; DelTx(P)) = star(p; De\(P w )). 



The equality of the Delaunay complexes now follows from Proposition 3.9 Equation 



The homeomorphism assertion follows from previous work |BG10| Theorem 5.1]. □ 



4 Algorithm 

In this section we introduce a Delaunay refinement algorithm which, while constructing a tangen- 
tial Delaunay complex, will transform the input sample set into one which meets the requirements 
of Theorem 1 3. 5 1 In particular we wish to construct a tangential Delaunay complex in which every 
m-simplex a is T -thick and for every p £ a, there is a 5 2 -power-protected Delaunay ball for a 
centred on T p A4. We demand 5 > (JoAoCj where e provides a strict upper bound on the radius of 
these Delaunay balls, and /ioe provides a lower bound on the shortest edge length of any simplex 
in DelTMi'P)- The constants <5o and fiQ are both positive and smaller than one. 

The algorithm is in the same vein as that of Boissonnat and Ghosh [BG10 , which is in turn 



an adaptation of the algorithm introduced by Li [Li03 . It is described in Section 4.2 after we 
introduce terminology and constructs which are used in the algorithm in Section |4.1| 

4.1 Components of the algorithm 

We now introduce the primary concepts that are used as building blocks of the algorithm. 



4.1.1 Elementary weight functions 

Elementary weight functions are a convenient device to facilitate the identification of simplices 
a that are not <5 2 -power-protected for 5 = SqL(o-). 

In order to emphasise that we are considering a function defined only on the set of vertices 
of a simplex, we denote by a the set {pq, . . . , p^} of vertices of a — [po ■ ■ ■ Pk\- We will call 
u a : a ~ > [0, oo ) an elementary weight function if it satisfies the following conditions: 
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1. There exists Pi & cr such that uj a (pi) £ [0, SoL(a)], and 

2. for all Pj <E ij\pi, u a (pj) = 0. 

For a given a = [p , ■ • • , Pfc] and elementary weight function co a , we define N(cr, cj ct ) as the 
set of solutions to the following system of k equations: 

11^ - -Pill 2 - 11^ - -Poll 2 = ^a{Pif -UJa(Pof- 

In direct analogy with the space N(a) of centres of a, the set N{a 1 u a ) is an affine space of 
dimension m — dimaff(cr) that is orthogonal to aff(cr). We denote by C(a,ui„) the unique point 
in N(a, uj a ) n afF(o-), and we define 

R(a,u a ) 2 = \\ Po - C(a,cj a )\\ 2 -u a (p ) 2 , 

where the notation is chosen to emphasise the close relationship with the circumcentre C (a) and 
circumradius R(a). The following lemma exposes some properties of R(a,u r7 ) in this spirit: 

Lemma 4.1 For a given a = [po, ■ ■ ■ ,Pk], with k > 1, and elementary weight function ov, we 
have: 

1. If a i < g then ui ai — u a \$ 1 is an elementary weight function, and 

J?(cti,w CTi ) < R(a,uj a ). 

2. A(a) < ^i?^,^). 

3. If T(cr) > 0, then 

< „, v <! + »?, 

i?(cr) 

With 7] = ^y. 

Proof 1. That cj cri is an elementary weight function follows from the observation that L(oi) > 
L(a). Since iV(CT, u ff ) C N(ai,co r7l ), the projection of C(«t, uj a ) into aff(oi) is C(<7i,u; CTl ). The 
result then follows from the Pythagorean theorem. 

2. Let e = [po>Pi] be the longest edge of a, and let c denote the projection of C(o-,cj a ) onto 
aff(e). Without loss of generality we assume that uj(p ) — 0. 

We have 

I bo - c|| 2 = ibi - c|| 2 - ^(pi) 2 

= ll(Pl -Po) - (c--Po)|| 2 - ^a(Pl) 2 

= A(a) 2 - 2(pi - po) • (c - Po) + Ibo - c|| 2 - oj a ( Pl ) 2 . 

Since p , Pi, an d c are colinear, we have 2(pi — p ) ■ (c — p ) = 2A(a) \\po — c||, and using the 
fact that Lo a (pi) < 8qL{&), we get 

lbo-c|| = ^^ 



> 



A{af 
(1 - 6 2 )A(a) 
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The result follows from the fact that R(a,uj a ) > ||po ~ c\\. 

3. Using the fact that oj a (p) = for all vertices p € a, except at most one, we get Pi € 
dB^k(C(a,ui a ), R(o~,uj a )) for all pi e a except at most one. 

Let t) = \\C(o~, uj a ) — C(<t)||, and assume, without loss of generality, that the vertex po € 
dB{C(o-,u a ),R{o-,u a )). Therefore, 

||C(<r)-p ||-»7 < II C(a,u a ) -poll < ||C(<r)-po||+»7 

R(a)+r) < R(cr,uj a ) < R(a) + r). (11) 

Since the point in N(a) that is closest to C(er, uj a ) is C(<r), and T(cr) > 0, we obtain the 
following bound using Lemma 4.1 from [BDG12J: 

5lL{af 



< 



2T(cr)A(cr) 



< Jr^R(cr), since L(a) < A(a) < 2R(a). (12) 



The result now follows from Eq. (11 1 and (12 1. □ 



If a = p * dp, and is an elementary weight function that vanishes on d p , then N(a, u a ) C 
N(a p ), but no point in N(a,uj r7 ) can be the centre of a <f 2 -power-protected Delaunay ball for a p 
for any 5 > <5o£(c). In other words, a and uj a define a quasi-cospherical configuration that is an 
obstruction to the power protection of o~ p at all points in N{a,uj a ). 

4.1.2 Quasicospherical configurations 

We now define the family of simplices that our algorithm must eliminate in order to ensure that 
the final point set has the desired protection properties. 

Recalling the definition Q of star(p), we have the following [ BG10I Lemma 2.7 (1)]: 

Lemma 4.2 Let V C A4 satisfy a sampling radius of e with respect to d R w such that e < 
rch(.A/f)/16. Then for all x e V k n (p) n T p M, we have \\p — x\\ < 4e. In particular, for all p e V, 
and every m-simplex a e star(p), we have R p (o~) < 4e. 



Since by Lemma 4.2 the Voronoi cell of p restricted to T p Ai is bounded, we get: 



Lemma 4.3 If e < rch 1 ^^ , then the combinatorial dimension of the maximal simplices in star(p) 
is at least m. 

We will always assume that V satisfies a sampling radius of e < rch 1 ^^ . If a is a maximal simplex 
in star(p), then Vrjv(ct) intersects T p M. at a single point. Indeed, since Vhn((t) C Vrn(p), by 
the convex set Vrn (ct) n T p M is bounded, and if it had a nonempty interior, then a 



4.2 



Lemma 

would not be maximal. Let a be a maximal simplex in star(p). Then, for all o~ m < er, the unique 
point in Vrjv(it) n T p M. . will be denoted by c p (a m ). We denote the radius of the circumscribing 
ball centred at c p (cr m ) by R p (a m ), i.e., R p {a m ) = \\p - c p (a m )\\. 

In our algorithm we will use the following complex, whose definition employs a particular 
elementary weight function: 

cosph^Op) = | a m+1 = Prn+1 * a m | a m £ star(p), R p {a m ) < e, 

a m is To-good, and 3 u a m+i with w a m+i \&m = 

and c p (a m ) e iV(a m+1 , ) I. (13) 
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The (m + l)-dimensional simplices in cosph 5o (p) are analogous to inconsistent configurations 
defined in (H£lOl |BG11J . 

Unless otherwise stated, whenever a m+1 = p m +\ * a m G cosph 50 (p) , with a m G star(p). the 



mention of w„m+i will refer to the elementary weight function identified in Equation (13). In 
particular, 

w CT m+i (pi) = for all pi G (J m+1 \ Pm+i, 

and 

^W^m+i) G [0,«5oi(a m+1 )] 

satisfies 

IMO -p|| 2 = ||cp(o- m ) -p m+ i|| 2 -awi(p m+1 ) 2 . 
We will exploit the following observations: 

Lemma 4.4 If a m+1 = p m +i * cr m G cosph 5 °(j>) with cr m G star(p), then 

i?(a m+1 ,^ m+1 ) <i?> m ) 

and 

Proof Since c p (cr m ) £ ./V(f7 Tn+1 ,av™+ 1 )j it follows that C(a m+1 , w CT m+i ) is the projection of 
Cp(a m ) into aff(cr m+1 ), and therefore R p (a m ) > R(a m+1 , u am +i ). The bound on A(cr m+1 ) now 
follows directly from Lemma |4.1| □ 



Boissonnat et al. |BG11| , using Lemma 4.2 showed that we can compute star(p) by computing 
a weighted Delaunay triangulation on T V M. of the points obtained by projecting V onto T p A4. 
Once star(p) has been computed, we can compute cosph 5 " (p) by a simple distance computation. 

The importance of cosph" 50 (p) lies in the observation that if an m-simplex o~ m G star(p) is not 
sufficiently power-protected, then there will be a simplex in cosph 50 (p) that is a witness to this. 
It is a direct consequence of the definitions, but we state it explicitly for reference: 

Lemma 4.5 If V is /2o£-sparse, and cosph* (p) — 0, then every a m G star(p) is <5 2 £tQ£ 2 -power 
protected on T p A4. 

4.1.3 Unfit configurations and the picking region 

The refinement algorithm, at each step, kills an unfit configuration by inserting a new point 
x = ijjp(x') where x' belongs to the so-called picking region of the unfit configuration, and ip p is 
the inverse projection defined in Equation pj. We use the term unfit configuration to distinguish 
the elements under consideration from other simplices. An unfit configuration <j) may be one of 
two types: 

Big configuration: An m-simplex cf> — <r rn in star(p) is a big configuration if R p {a m ) > e. 

Bad configuration: A simplex <f> is a bad configuration if it is To-bad and it is either an m- 
simplcx <fi = a m G star(p) that is not a big configuration, or it is an (m + l)-simplex 
(j) = a m+1 G cosph 50 (p). 



5.2 



Lemma 



5.13 



that in fact every (m -I- l)-simplex in cosph 5 °(p) is a 



We will show in Section 
bad configuration. 

The size of the picking region is governed by a positive parameter a < 1 called the picking 
ratio. 
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Definition 4.6 (Picking region) The picking region of a bad configuration, a m € star(p) or 
Pm+i * °~ m G cosph s °(p) with a rn € star(p). denoted by P(<r m ,p) and P(cr m+1 ,p) respectively, is 
defined to be the m-dimensional ball 

B RN (c p (a m ),aR p (a m ))n T p M. 

We choose a point in the picking region so as to minimize the introduction of new unfit 
configurations. We are able to avoid creating new bad configurations provided that the radius 
of the potential configuration is not too large. To this end, we introduce the parameter f3 > 1. 

Definition 4.7 (Hitting sets and good points) Let <f> — a m G star(p) or <fi = q * a m g 

cosph 5 °(p) with a m G star(p), and x = ip p (y) where y £ P(0,p). A set <r C f of size fc, with 
k < m + 1, is called a hitting set of a; if 

a. t = i * cr is a fc-dimensional To-flake 

and there exists an elementary weight function uj t satisfying the following condition: 

b. R(t,u t ) < /3R p (a m ) 

The elementary weight function uj t is called a hitting map, and we sometimes say a hits x. 
A point x — ip p (y), where y € P(4>,p), is said to be a good point if it is not hit by any set 
crcP with \a\ < m+ 1. 

A simplex cr which defines a hitting set of x, is necessarily To-good. This follows from the 
requirement that x * a be a To-flake. 



4.2 The refinement algorithm 

In this section, we show that we can refine an e-net of M. so that the simplices of the Delaunay 
tangential complex of the refined sample DelrAi('P) are power-protected. An e-net is a point 
sample P cM that is an e-sparse e-sample set of M. for the metric o\n . One can obtain an 
e-net by using a farthest point strategy to select a subset of a sufficiently dense sample set. We 
will assume that we know the dimension m of the submanifold A4 and the tangent space T p A4 
at any point p in M. 

The algorithm takes as input Po, an e-net of M, and the positive input parameters e, Tq, 
a < |, P > 1 and <5o < 4. The algorithm refines the input point sample such that: 

(1) The output sample P D P is an /io e_s P arse e-sample set of with respect to g? k « , where 
Mo - g. 

(2) For all p £ V, every m-simplex a m 6 star(p; DelTvw(P)), f m is To-good and (5Q/io e2 -P ower 
protected on T p M. 

The algorithm, described in Algorithm [T| applies two rules with a priority order: Rule (2) is 
applied only if Rule (1) cannot be applied. The algorithm ends when no rule applies any more. 
Each rule inserts a new point to kill an unfit configuration: either a big configuration or a bad 
configuration. 

A crucial procedure, that selects the location of the point to be inserted, is Pick_valid, given 
in Algorithm[2] Pick_valid(0,p) returns a good point ip p (y) where y G P((j),p). 

The refinement algorithm will also use the procedure Insert (p), given in Algorithm p] 



5 Analysis of the algorithm 



We now turn to the demonstration of the correctness of Algorithm [T] In Section 5T we show 
that the algorithm must terminate, and in Section [5. 2| we show that the output of the algorithm 
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Algorithm 1 Refinement algorithm 

Input e-net Vo of M, and input parameters r , a and So; 
Initalize V Vo, and calculate DbItmCP); 
Rule (1) Big configuration (e-big radius): 

ifBpeV such that 3 a m e star(p) with R p (a m ) > e, 

then Insert(^ p (cj,((T m ))); 
Rule (2) Bad configuration (To-bad): 

if 3 p e V and 3 a rn e star(p) s.t. <7 m is r -bad, 

then Insert(Pick_valid(er m ,p)); 

if Bp e V and 3 a m+1 € cosph 5 °(p) s.t. cr m+1 is r -bad, 
then Insert(Pick_valid(cr m+1 ,p)); 
Output DelTx('P) = Upg-p star(p); 



Algorithm 2 Pick_valid(<7, p) 

1 1 Assume that a is either equal to a m € star(p) 

//or cr m+1 = p m+1 * cr m e cosph 5o (f>) with cr m G star(p) 

Step 1. Pick randomly y <= P{<r m ,p) (or P(<r m+1 , p)); 

/ / Recall that i/> p projects points from TpA^ onto .M along iVpAl 

Step 2. x «- V P (y); 

Step 3. Avoid hitting sets: 

1 1 \o\ denotes the cardinality of a 

if 3 <7 C V, with | cr | < m + 1, which is a hitting set of x, 
then discard x, and go back to Step 1; 
Step 4. Return x; 



Algorithm 3 Insert(p) 
Step 1. Add p to V; 

Step 2. Compute star(p) and cosph 5o (p); 

Step 3. For all x e "P \ {p}, update star(x) and cosph 5o (x); 
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meets the requirements of Theorem |3.5[ In order to complete the demonstrations we impose a 
number of requirements on the input parameters, listed as Hypotheses HO to H5 below. 

Recall that our input parameters are the following positive numbers: e, which is the sampling 
radius and sparsity bound satisfied by Vq, the input e-nct sample set; So, which is used to describe 
the amount of power-protection enjoyed by the m-simplices in the final complex; Tq, which is 
used to quantify the quality of the output simplices; /3, which is used to describe an upper bound 
on the radius of the bad configurations that we will avoid; and a, which governs the relative size 
of the picking region. 

It is often convenient to represent the sampling radius by a dimension-free parameter that 
has the reach of the manifold factored out. We define 

e 

e = 



TCh(M)' 

The volume of the m-dimcnsional Euclidean unit-ball is denoted V m . In order to state the 
hypotheses on the input parameters, we use some additional symbols: 

1 



2 4 (2 4 + l) 2 ' 
J B = 4 + 2(l + 2 7 3 2 /? 2 ) 2 , 



1-2% = 



as well as £, E, and D. The term £ is introduced in Lemma 5.5 and depends on m and rch(.M) 



and the term E, defined in Equation (17 1, depends on £ and (3. The symbol D is introduced in 
Lemma [5. 8| where it is said to depend on m and (3. 

In order to guarantee termination, we demand the following hypotheses on the input param- 
eters: 



HO. a < 1/2 

m - 13 ~ (l-5g)(l-a-4.5e ) 



m. si < ix +1 



HA. e < min 



m \2(/3+/3')rch(VW)' 8/3 / 



To meet the quality requirements of Theorem |3.5| we demand an additional constraint on the 
sampling radius: 

H5. e < ^Lt" 



l.lxlO 9 

The make use of the following observation: 
Lemma 5.1 From hypotheses HQ to HA we have e < eo and <5 2 < 2 4 eo, and 

(l-<Sg)(l-a-4.5e )e e de f . n 
1 > 9 - Moe- (14) 
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Proof From HI we have (3 > 2 and using the fact that B > (3 4 and H2 we have Tq < pVj- 
And using the fact, from "H4, that 



- r o' i+1 r o 1 

e < — < — < = e n . 

" 8/3 ~ 8(3 2 4 (2 4 + l) 2 



Similarly the bound on 5q follows from H'S. 



Inequality (14 1 follows from HO and the definition of Cq. 



□ 



From Equation (14) we can see that we require (3 > 4.5. Given a satisfying HO, and a valid 



choice for (3, the hypotheses H2 to H4 sequentially yield upper bounds on the parameters Tq, 
5q, and e; we are able to choose parameters that satisfy all of the hypotheses. 
The main result of this section can now be summarised: 



Theorem 5.2 (Algorithm guarantee) If the input parameters satisfy hypotheses HO to H5, 
then Algorithm [l] terminates after producing an intrinsic Delaunay complex T)e\j^(V) that tri- 
angulates M.. 

5.1 Termination of the algorithm 

This subsection is devoted to the proof of the following theorem: 

Theorem 5.3 (Algorithm termination) Under hypotheses HO to HA, the application of Rule (1) 
or Rule (2) on a big or a bad configuration <f> always leaves the interpoint distance greater than 



and if is a bad configuration then there exists x € P{<p,p) such that ip p (x) is a good point. 
Since A4 is a compact manifold this implies that the refinement algorithm terminates and returns 
a point sample V which is an /2o£-sparse e-sample of the manifold M.. 

We will prove that at every step the algorithm maintains the following two invariants: 

Sparsity: Whenever a refinement rule inserts a new point x = if> p (y), the distance between x 
and the existing point set V is greater than jX e. 

Good points: For a bad configuration <j> refined by Rule (2) , there exists a set of positive volume 
G C P(cj),p) such that if x € G, then tp p (x) is a good point. 

The Termination Theorem |5.3| is a direct consequence of these two algorithmic invariants. We 
first prove the sparsity invariant in Section |5.1.1| using an induction argument that relies on the 
fact that the algorithm only inserts good points. The existence of good points is then established 
in Section |5.1.2| using the sparsity invariant and a volumetric argument. Termination must follow 
since M is compact and therefore can only support a finite number of sample points satisfying 
a minimum interpoint distance. 



5.1.1 The sparsity invariant 

The proof of the sparsity invariant employs the following observation, which serves to bound the 
distance between a point inserted by Rule (2) and the existing point set: 
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Lemma 5.4 Assume Hypotheses HO to HA. Let <f> = a m £ star(p) or <fi — p m +i * u m £ 
cos]A\ So (p) be a bad configuration being refined by Rule (2). Then for all x £ P(<j),p) we have 

d R ,(c p (a m ),^ p (x)) < (a + A.5i )R P (a m ) 

and 

d R *ty p (aO,7>) > (1 - a - 4.5i )R p (a m ) > ~ J ^~ A - 
Proof Using the facts that a < |, and e < eo, and R p (a m ) < e, we have that for all a; £ P(<fi,p) 

\\p- x \\ <{ l + a)R p (a m )<^<^<\, 

and so we may apply Lemma |B.2| to get 

. , ... ^ 2||p-:r|| 2 2(l + a) 2 i? p (a m ) 2 ^ , c . D , m , 



and 



\\c p {a m ) - *P p {x)\\ < \\c p (a m ) - x\\ + \\x - ijj p (x)\\ 

< (a + 4.5e ) R P (o- m ). (15) 



Let S p = dB R N(c p (a m );R p (a m )). From Eq. g5j we have for i e P(^,p) 

> (l-a-4.5e ) R p (cr m ) 
R p (a m ) 
3 ' 

where the final inequality follows from HO and the definition of eo- D 

We introduce some additional terminology to facilitate the demonstration of the sparsity 
invariant. An abstract simplex in the initial sample set a c Po is called an original simplex, 
otherwise a C P is called a created simplex. 

Let </> be an unfit configuration that was refined by inserting a point x. We say that x created 
a if x £ a and x is the last inserted vertex of the simplex a, i.e., a \ {x} already existed just 
before the refinement of the unfit configuration <f>. The unfit configuration (f> is called the parent 
of a and will be denoted p(cr). 

Let a denote the simplex being refined by the refinement algorithm. We will denote by e(cr) 
the distance between the point newly inserted to refine a and the current sample set. 

The sparsity invariant is demonstrated by induction. We use a case analysis according to 
the type of unfit configuration being refined; it is necessary to consider sub-cases. The in- 
duction hypothesis is employed only in the sub-case Case 2(b) (ii) and the implicit similar 
Case 3(b) (ii). The base for the induction hypothesis, i.e., the insertion of the first point, 
cannot involve Case 2(b) or Case 3(b). 

Case 1. Let <j) = a m £ star(p) be a big configuration being refined by Rule (1). 



4.2 



Since Po (C P) is an e-net, we have from the fact that e < eo < jq and Lemma 
R p (a m ) < 4e. Rule (1) will refine a by inserting ip p {c p (a m )) . Using the fact that e < eo < 
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Yg, R p {a m ) < 4e and R p (a m ) > e (since a m is being refined by Rule (1)), and Lemma B.2 
the distance between i(j p (c p (a m )) and any vertex inserted before ip p (c p (a m )) is not less than 

R p (a) - \\c p (a) - Mcp(*))\\ > W ~ 

> (l-8e )e 
e 

>2< 



which establishes the sparsity invariant for this case. 



Case 2. Consider now the case where <f> = p m +i * o m € cosph 50 ^), with er m € star(p), is being 



refined by Rule (2). In this case, recalling Lemma 2.5 we have 



• R p {(T m ) < e, and 

• there exists a face of <p that is a r -flake. 

Let <Ji C <fi denote a face of <j> that is a r -flake. We have to now consider two cases: 

(a) <j\ is an original simplex 

(b) a\ is a created simplex 

Case 2(a). If o\ is an original simplex then a\ C Po, an d since is an e-net, L(ai) > e. Since 
a flake must have at least three vertices, cri and u m must share at least two vertices, and 
therefore R(a m ) > e/2. 



Let cc = i/j p (x') be point inserted to refine where x' € P(<p,p). Using Lemma 5.4 and the 
fact that R(a m ) > e/2, we therefore have 

d RN (x, V)>(l-a- 4.5 e )R P {o- m ) 
> {l-a-4.5i )R(a m ) 
(l-a-4.5e )e 
2 



where the final inequality follows from Inequality ( 14 1 . Hence the sparsity invariant is 



maintained on the refinement of if o\ is an original simplex. 

Case 2(b) We will now consider the case when g\ is a created simplex. We denote by p(ui) 
the parent simplex whose refinement gave birth to o\ . 

We will bound the distance between x = ip v (x'), where x' € P{4>,p), and the point set V '. 
Let x* denote the point whose insertion killed p(cti). By definition x* is a vertex of a±, 
and hence also of <f> since o\ < <f>. We distinguish the following two cases: 

Case 2(b) (i) Suppose p{a{) was a big configuration refined by the application of Rule (1). 
According to Case 1, the lengths of the edges incident to x* in u\ are greater than e/2. 
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Therefore 

d RN (x,V) > (1 - a- A.5e )R p (a m ) by Lemma [Q| 

(1-$)(1 -a-4.5e )AU) , __ 

> ^ ^ u; by Lemma [44] 

> ( 1 -^( 1 "7 4 - 5g °^ since crj < 

(l-^)(l-a-4. 5g0 ) e Cage 

4 J 

> fl^e Inequality [T4j 

and the sparsity invariant is maintained. 

Case 2(b) (ii) Suppose p(<7i) was a bad configuration refined by Rule (2). Thus p(cti) was 
either an m-simplex er™ €E star(q) or an (m + l)-simplex q m +i * cr™ G cosph 5 °(g) with 
a? e star(g). 



Consider the elementary weight function uj ai — w^lcr'i , where is the weight function ( 13 1 
identifying </> as a member of cosph* 50 ^). From Lemma 4.1 1), and Lemma 



4.4 



we have 

that R p (a m ) > R(ai,uj ai ). We also have that R(ai,oj ai ) > f3R q (a?). Indeed, otherwise 
°i \ i x *} would be a hitting set for x* , contradicting the hypothesis that p(<7i) was refined 
according to Rule (2) by the insertion of a good point x* . Thus we have 

d RN (x, V) > (1 - a - 4.5e )i? p (cr m ) 

> (1 -a-4.5e )i?((Ti,a; (Tl ) 

> (l-a-4.5e )/3i?> 2 m ) 

. (l-<5 2 )(l-a-4.5e )/3A(a 2 m ) „ 

> — — — x -^-L from Lemma [4~4l 



2 

> A(cr™) from Hypotheses HI on /3 

> Aoe, 

where the last inequality follows from the induction hypothesis. Again the sparsity invariant 
is maintained after refinement of <p. 

Case 3 The proof for the case of a bad configuration <\> = a ,n G star(p) to be refined by Rule (2) 
is similar to Case 2, and the lower bound on the interpoint distances is the same. 

This completes the demonstration of the sparsity invariant. 
5.1.2 The good points invariant 

We will now show that the good point invariant is maintained if <fi is a bad configuration being 
refined by Rule (2). Without loss of generality, we will assume that <f> is either equal to a rn G 
star(p) or to q * a m G cosph 5 °(p), with a rn G star(p). 



Recall the picking region P(cj>,p) introduced in Definition 4.6 We will show that there exists 
y G P(4>,p) such that x = ip p (y) is a good point. Let Y C P((j>,p) be the set of points that ijjp 
maps to a point with a hitting set: 

Y = {y G P((f),p) | ip p (y) is not a good point}. 
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We will show that the volume of P(4>,p) exceeds the volume of Y. To this end, we will first 
bound the number of simplices that could hit some point in ip p (Y). Then we will bound the 
volume that each potential hitting set can contribute to Y. 

In order to bound the number of hitting sets, we will use the sparsity invariant together with 
the following lemma [BG10, Lemma 4.7] to bound the number of points that can be a vertex of 
a hitting set: 

Lemma 5.5 (Bound on sparse points) For a point peM and R > 0, let V be a maximal 
set of points in B r n\ m (p, R) such that the smallest interpoint distance is not less than 2r. There 
exists £ that depends on m and rch(.M), and A that depends on m, such that if R + r < £, then 

We obtain the following bound on the number of hitting sets: 

Lemma 5.6 Let S(<fi) denote the set of simplices contained in B + n V that can hit a point in 
ip p (Y). Then 

zpm+l 



where 



E= 2 



1 + A£ 



(18(a + 2/3'+6.5e ) + l) r 



(17) 



Proof Suppose a C V is a hitting set of a point x = ip p (y), where y € P((f>,p), with \a\ — 
k and k < m + 1. Let a = x * a, and let u>s- denote the corresponding hitting map (see 
Definition 4.7). Therefore, we have R{a,uj ^) < (3R p (a m ), and it follows from Lemma 4.1 ^2) that 

A 0) < T3§^-R P (0- Thus from Lemma 
B R N(c p (o- m ),r~), where 



5.4 



del 



4.5e R p {a m ) + aR p (a m ) 



and the Triangle inequality we have a C B = 

2/9 



R P (a m ). 



Let c = ipp{c p {a m )). Then using Lemma B.2 and the fact that i? p (cr m ) < e < e rch(7W) we 
have 

Using (5q < 2 4 eo from Lemma 5.1 and j3' = 1 _^ 4 - , and i? p (er m ) < e, we find 

IMO - c|| + r- < aiZp(a ro ) + 6.5e R P (a m ) + 

< (a + 2^' + 6.5e )e 
= f ii. 



Thus B C B + = £ B^n (c, R), and y € Y if and only if there exists cr c B + n "P such that a hits 
^p(2/)- 
Us 

and observe that 



Using Lemma 5.5 we will bound the number of sample points in B + n V . Set r = ^ 



18 



R 



18 



a + 2/3' + 6.5 e U < (2/3' + l)e < f , 
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by Hypothesis T-LA. The sparsity invariant and Lemma [5 . 5 1 then yields 

|B+np|s ^ x( «^™ + l y^ f 

Since the number of fc-simplices is less than ) k+1 , and the maximum dimension of a hitting 
set is m, we have |S((/>)| < E ^ . □ 

We now turn to the problem of bounding the volume of Y. We will consider the contribution 
of each a G S (</>). The following definition characterises the set of points in M. that can be hit 
by o-: 

Definition 5.7 (Forbidden region) For a /c-simplex a with vertices in M. with k < m and 
parameter t < e, the forbidden region, F(a,t), is the set of points x £ Ai such that ct\ — x * a 
satisfies the following conditions: 

• L(ax) > § 

• (jj is a r -flake 

• there exists an elementary weight function u> ai s.t. R(o~i, uj ai ) < fit 

We will use the following lemma, which is proved in Appendix [C] It bounds the volume of the 
set of points that can be hit by a given simplex: 

Lemma 5.8 (Volume of forbidden region) Let a be a fc-simplex with vertices on M. and 
k < m. If 

2 - 1 ^ min { 4/ 3rch(^) ' and 

3. <5 2 <min{r™ +1 ,i}, 
then 

vol(F(a,t)) <Dr R(o-) m 7 

where D depends on m and fi. 



Lemma |5.8| together with Lcmma |5.6[ yields a bound on the set of points Y in the picking region 
that do not map to a good point : 

Lemma 5.9 The volume of the set Y C P(<j),p) of points that do not map to a good point is 
bounded as follows: 

vol(F) < E m+1 fi m DT R p (a m ) m . 
Proof Let to = R P {c m ) < £■ For a given a E S(</>), let Y a C Y be the set of points y for which 



a hits x = ip p (y). Then from Hypotheses HO to "H4 and Lemma 5.8 we have 



vol(F CT ) <vol(TT p (F(a,t ))) 

< vo\(F(a,t )) since tt p is a projection map on T p M 

< DT R(a) m . (18) 
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Let a i — x * <7, and let uj ai be the corresponding hitting map. From the definition of hitting 
sets and hitting maps, we have R p (a m ) < e, and R{a\ , uj ai ) < (3R p (a m ) and a is Tp-thick. Define 



we have 



Then, using Lemma 4.1 (3) and the fact that R(a,uj a ) < R(ai,uj ai ) < f3R p (a r ' 



< R(o-,oj a ) (l 

< 2R{a.u a ) 

< 2/3R p (a m ). 



T(a) 
S 2 o 



since T(cr) > Tg > 



§2 1 

since — < T < - from Hyp. H2, U3 



(19) 



The inequalities ( 18 1 and ( 19 ) together yield 



vol(r CT ) < 2 m l3 m DT () R p (a r 
and so using Lemma [5. 6| we have 



vol(y)=vol |J Y a 

< £ voi(y CT ) 

cr£S(0) 

< E m+1 (3 m DT R p (cr m y 



By the definition of the picking region, we have that 



□ 



vol(P(0,p)) = K l a m i?p(^ m ) m . 
By Hypothesis H2, E m+1 /3 m DT Q R p {a m ) m is less than vo\(P((f>,p)), the volume of the picking 



region of <fi. Thus with Lemma 5.9 this proves the existence of points y in the picking region 
P(4>,p) of 4> such that ippiu) is a good point. 
The proof of Theorem |5.3| is complete. 



5.2 Output quality 

We will now show that if Hypothesis H5 is satisfied, in addition to Hypotheses HO to RA, then 
the output to the refinement algorithm will meet the demands imposed by Theorem |3.5| thus 
yielding Theorem |5.2| 

The main task is to ensure that every m-simplex in T)g\tm ("P) has, for each vertex, a SQ^e 2 - 
power-protected Delaunay ball centred on the tangent space of that vertex. This is achieved in 
two steps. First we establish conditions to ensure that cosph* 50 ^) = for every p £ V . As noted 



by Lemma 4.5 this ensures that every simplex in star(p) has a ^Q/2ge 2 -power-protected Delaunay 
ball centred on T p A4. Next we show conditions such that if a rn £ DelT_vi (T 5 ), then a m E star(p) 
for every vertex p € o~ m . In each step the required conditions impose an additional constraint on 
the sampling radius, and this leads to Hypothesis H5. 

As a starting point, we observe the following direct consequence of the Termination Theo- 
rem 15.31 
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Corollary 5.10 Under Hypotheses HO to %4, for all p £ V, the output of the algorithm satisfies 
the following: 

1. er m 6 star(p) => R p (a rn ) < e and cr m is a To-good simplex, and 

2. all cr m+1 € cosph 5o (p) are r -good. 

We will show that for an appropriate sampling radius, there cannot be a To-good simplex in 
cosph °{p). We exploit the following bound on the thickness of a small (to + l)-simplex: 

Lemma 5.11 (Small (m + l)-simplices are not thick) Let a m+1 be an (to + l)-simplex 
with vertices in M and A(cr m+1 ) < rch(A^). For distinct vertices p, q G <j m+1 define 9 = 
Z(&S(a q ),T p M). Then 

'A(<7 m+1 ) 



T(a m+i ) < 



sin 9 



v 2rch(X) 

Proof We will bound the altitude -D(g, u m+1 ). Let I be the line through p and g. Using 
Lemma B.l and the fact that Z(aff (cr q ), T p M) = 9, we get 

D(q,<T m+1 )=d RN (q, a $(<T q )) 

= sinZ(£, aff(<r g )) xd R N(p 7 q) 

< (sin Z(£,T p M) + sin Z[afi(a q ),T p M)) x d R «(p,q) 
' d R N (p, q) 



< 



< 



2rch(7W) 

A((7 m+1 ) 

2rch(M) 



sin J x d R « (p, q) 
+ sin6M x A(d m+1 ). 



Therefore we have 



T(o- m+i ) < 



A(cr m+1 ) 

2rch(X) 



sm ( 



□ 



Also, Whitney's Lemma 2.1 implies that a T -good simplex in star(p) makes a small angle with 
the tangent space at p: 

Lemma 5.12 If a m e star(p) is T -good with R p (a m ) < e, then 

s[n6< f^iMY 

where 6 = Z(aS(a m ),T p M). 

Proof Let £ = max ief , m g^n (x, T p M) where a; is a vertex of er m . From Lemma 



B.l 



we have 



£ = max d R n(x,T p M) 



d R «(P, x) 2 
2rch(.M) 

A(q m ) 2 
2rch(7W)' 



< max 



< 
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Using Lemma [2TT] and the facts that R{a m ) < R p {a m ) < e and T(a m ) > T^ (since a m is a 
r -good simplex), we have 



sin 9 < 



< 



2( 



< 



T(cr m )A((T m ) 

A{a m ) 
T(a m )rch{M) 
2e 

T™ rch(X) 



since A(a m ) < 2i?(a m ) < 2e. 



□ 



Using Lemmas 



5.11 



and 



5.12 



we get that no (m + l)-dimensional simplices in cosph °(p) can 



be To-good when e is sufficiently small: 

Lemma 5.13 (cosph 50 (p) simplices are To- bad) Let a m+l = p m+1 * a rn £ cosph" 50 (p) with 
a m £ star(p). If 



e < 



^2m+l 




and 5g < \, then T(a m+1 ) < T^ 



m+l 



Proof By Lemma |4.4| 



since R p (a m ) < e and (S§ < ~. 



A(a 



m+l\ 



< 



,-i? p (a m ) < 4e, 



Then from Lemmas |5.11 and |5.12| we get 

T(a m+1 ) < 



< 



< 



/A(a m+1 ) 
V2rch(A4) 

2e 
rch(A^) 
4e 



1 



r "rch(X) 



- sin 
1 

1 

since r < 1 



<r, 



m+l 
) 



from the hypothesis on e. 

We emphasise the consequence of Lemma |5.13| 



□ 



Corollary 5.14 If 8% < \ and 



e < 







and all the simplices in cosph 5 " (p) are To-good, then 

cosph 50 (p) = 0. 

Now we proceed to the second step of the analysis. Assuming that cosph 50 (p) = for all p 
in V, the following lemma says that if a £ star(p), then also a £ star(g) for every vertex q £ a, 
provided the appropriate constraints are met. 
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r2 .-.2pm 



Lemma 5.15 Let V be a /ioe-sparse e-sample of .M with /io < 1 independent of e. We further 
assume 5 a < 1 and 

(1) for all p € "P, every er m G star(p) is a To-good simplex with R p (a m ) < e, and 

(2) for all cosph s °(p) = 0. 

If 

36 ' 

then star(p) = star(p; Delrx (V)) for all p in P. 

Proof For p G "P, let cr m € star(p) and q p) be a vertex of a m . We will show that a m is 
also in star(q). 

Let 9 = maxZ(aff (<7 m ), T x Ai) where the max is taken over the vertices x of a m . Since 

" 36 4 ' 

Lemma |5.12| yields 

. 2e def _ 1 

smfl < = cie < -. 

- pm 2 

It follows that cos# > V3/2 and so 

tan# < 2cie. 

Recall that N(a m ) denotes the affine space orthogonal to aff (a m ) and passing through C(er m ). 
Let c be the unique point in N(a m ) H T q A4, and let R — d R N(c,p). 
Using the fact that Z(aff (cr m ), T q M) < 6, we have 

d RN (C(a m ),c) < R(a m ) t&n6 < 2 Cl eR(a m ), 

and likewise 

d R N(C(a m ),c p (a m )) < 2 Cl ~eR(a m ). 

It follows that R< {1 + 2c 1 i)R(a m ), and d RN (c p (a m ), c) < A Cl eR(cr m ). From the above obser- 
vations, and using the fact that R(a m ) < R p (a m ) < e, we get 

B R «(c,R) C B RN (c p (a m ),(l + 6 Cl i)R(a m )) 
C B R „(c p (a m ),i? p ( ( T m ) + 6c 1 ee). 

Since cosph 5 °(p) = 0, and V is /ioe-sparse, we have that a m is <5o/ioe 2 -power protected on 
T p A4 (Lemma 4.5). This means that 

B RN {c p {a m ), R p (a m ) + A) n (V \ a m ) = 0, 

where 

A = ^ R p (a^Y + 8l~pl^ ~ R p {a m ) 

> goAje 

> „ — C2£. 
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Since Qc\te < C2£, by our hypothesis on e, we have 

B RN (c,R) C B RN (c p (a m ),R p (a m ) + A), 
and thus the m-simplex er m belongs to star(g). 



□ 



The consequence of Lemma 5.15 together with Lemma 4.5 is that every m-simplex in DeLrx('P) 
has, for each vertex, a J^/ipe^-power-protected Delaunay ball centred on the tangent space of that 
vertex: 

Corollary 5.16 Let V be a /2 e-sparse e-sample of M. with jl Q being independent of e. Under 

in star(p;DelrAi(7 3 )) are 
€ star(p; DelTAi(7')) there exists a c p (<7 m ) € 



the hypotheses in Lemma 5.15 for all p 6 V, all the m-simplices a r 
#o/ioe 2 -power protected on T p M. I.e, for all a r ' 
N(a m ) n T p X such that for all q e V \ a m 



d R N{q,c p {a m )) 2 > d M N(p,c p (a m )) 2 +5fa 2 e 2 . 

We are now in a position to show that Hypothesis "H5, when added to Hypotheses HO to H4, 
results in the output of the algorithm meeting the demands of Theorem |3.5| 

Recalling that /zq = |, Hypotheses W3 yields the following consequence of T-L5: 



X2 j--2ra 
°0 1 



e < 



1.1 x 10 9 



< min 



F 2m+1 
1 



3G 



1.5 x 10 6 



In other words, the sampling radius bounds demanded by Corollary |5.14[ Lemma 5.15| and The- 
orem |3.5l are all simultaneously satisfied. Corollary |5. 1 0| together with Corollary |5. 14 ensure that 
the hypotheses of Lemma 5.15 are satisfied, and so it follows that the m-simplices of DelxA-f (T'O 



are power-protected as described by Corollary |5.16| Thus all the requirements of Theorem |3.5 
are satisfied, and we obtain Theorem 15.2 



6 Conclusions 

We have described an algorithm which meshes a manifold according to extrinsic sampling condi- 
tions which guarantee that the intrinsic Delaunay complex coincides with the restricted Delaunay 
complex, and that it is homeomorphic to the manifold. The algorithm constructs the tangential 
Delaunay complex, which is also shown to be equal to the intrinsic Delaunay complex, and in 
this way we are able to exploit existing structural results |BG11| to obtain the homeomorphism 
guarantee. 

This approach relies on an embedding of Ai in M. N . In future work we aim to develop algo- 
rithms and structural results which enable the construction of an intrinsic Delaunay triangulation 
in the absence of an embedding in Euclidean space. 
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A An obstruction to intrinsic Delaunay triangulations 

When meshing Riemannian manifolds of dimension 3 and higher using Delaunay techniques, 
flake simplices pose problems which cannot be escaped simply by increasing the sampling den- 
sity. In particular, developing an example on a 3-manifold presented by Cheng et al. [CDR05 , 
Boissonnat et al. [BGO09I Lemma 3.1] show that the restricted Delaunay triangulation need not 
be homeomorphic to the original manifold, even with dense well separated sampling. 

In this appendix we develop this example from the perspective of the intrinsic metric of the 
manifold. It can be argued that this is an easier way to visualize the problem, since we confine our 
viewpoint to a three dimensional space and perturb the metric, without referring to deformations 
into a fourth ambient dimension. This viewpoint also provides an explicit counterexample to the 
results announced by Leibon and Letschcr [LLOO : In general the nerve of the intrinsic Voronoi 
diagram is not homeomorphic to the manifold. The density of the sample points alone cannot 
guarantee the existence of a Delaunay triangulation. 

We explicitly show how density assumptions based upon the strong convexity radius cannot 
escape the problem. The configuration considered here may be recognised as essentially the same 
as that which was described qualitatively in Section |2.4.3| but here we consider the Voronoi 
diagram rather than Delaunay balls. We work exclusively on a three dimensional domain, and 
we are not concerned with "boundary conditions"; we are looking at a coordinate patch on a 
densely sampled compact 3-manifold. 

A.l Sampling density alone is insufficient 

We will now construct a more explicit example to demonstrate that the problem of near- 
degenerate configurations cannot be escaped with the kind of sampling criteria proposed by 
Leibon and Letscher [LL00_. 

Leibon and Letscher |LL0Q[ p. 343] explicitly assume that the points are generic which they 
state as 

Definition A.l The set V C A4, is generic if A4 is an m-manifold and m + 2 points never lie 
on the boundary of a round ball. 

Here a round ball refers to a geodesic ball. This definition of genericity is natural, and corresponds 
to Delaunay's original definition |Del34 , except Delaunay only imposed the constraint on empty 
balls. A question that Delaunay addressed explicitly, but which was not addressed by Leibon 
and Letscher, is whether or not such an assumption is a reasonable one to make. Delaunay 
showed that any (finite or periodic) point set in Euclidean space can be made generic through an 
arbitrarily small affine perturbation. That a similar construction of a perturbation can be made 
for points on a compact Riemannian manifold has not been explicitly demonstrated. However, in 
light of the construction we now present, it seems that the question is moot when m > 2, because 
an arbitrarily small perturbation from degeneracy will not be sufficient to ensure a triangulation. 

Leibon and Letscher proposed adaptive density requirements based upon the strong convexity 
radius. These requirements arc somewhat complicated, but they will be satisfied if a simple 
constant sampling density requirement is satisfied. Exploiting a theorem [Cha06, Thm. IX.6.1], 
that relates the strong convexity radius to the injectivity radius, inj(Al), and a positive bound 
on the sectional curvatures, they arrive at the following: 

Claim A. 2 Q LLOOl Lemma 3.3]) Suppose /Co is a positive upper bound on the sectional cur- 
vatures of A4, and 




(20) 
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If V is an ?7(A4)-sample set for M. with respect to d,M> then |Deljvi("P)| = A4. 

In fact, we will show that no sampling conditions based on density alone will be sufficient to 
guarantee a homeomorphic Delaunay complex in general, even when a sparsity assumption is 
also demanded. An i-net is an e-sparse, e-sample set. We will show: 



Theorem A. 3 With rj(A4) as defined in Equation (20 1, for any e > 0, there exists a compact 
Riemannian manifold M, and and a finite set V C A4, such that V is an (e?7(.A/f))-net for A4, 
with respect to the metric d^vi, but Del»('P) is not homeomorphic to M. 



A. 1.1 A counter-example 

We will construct the counter-example by considering a perturbation of a Euclidean metric. This 
is a local operation, and the global properties of the manifold are only relevant in so far as they 



affect ri(Ai) of Equation (20). We may assume, for example, that the manifold is a 3-dimensional 
torus M. = S 1 x S 1 x S 1 , initially with a flat metric. 

Thus assume there is some eo such that any compact Riemannian manifold may be triangu- 
lated by the intrinsic Delaunay complex when V is an eo?y(A / l)-net. For convenience, we choose 
a system of units so that eo?/(.M) = 1. We will first construct a point configuration and metric 
perturbation that leads to a problem, and then we will show that the sampling assumptions are 
indeed met. 

We introduce a number of parameters which we will ma- 
nipulate to produce the counter-example. We are exploiting 
the fact that the genericity assumption allows configurations 
that are arbitrarily close to being degenerate. The assumed 
eo has been fixed. 

We will work within a coordinate chart on M. , where the 
metric is Euclidean. We will perturb this metric by con- 
structing a metric tensor g, and we will denote by M. the 
manifold with with this new metric. 

Consider points u y v,w,p in the xz-plane arranged with 
u and v at ±a on the z axis, and w and p at ±(a + £) on the 
x axis, with a = § , and < £ < ro7, where tq and 7 will 
be specified below. The Voronoi diagram of these points in 
the xz-plane is shown in Figure [2] The main point here is 
that the Voronoi boundary between Vm ( u ) an d Vm ( v ) ma Y 
be arbitrarily small with respect to the distance between the 
sites, i.e., £ will be very very small. 

The three dimensional Voronoi diagram is the extension 
of this in the horizontal y-direction, so that every cross- 
section looks the same. Note that since the points are not co-circular, they do not represent 
a degeneracy by Delaunay 's criteria |Del34| , but this is irrelevant; we will also argue that the 
points will not represent a degenerate configuration with respect to the new metric. 

We now introduce a small localized metric perturbation so as to change the Voronoi diagram 
near the origin. For example, we can demand that the matrix of the metric tensor in our 
coordinate system has the form 





u 






• 




V 




• w 




• 






V 





Figure 2: A vertical slice: the xz- 
plane of the initial Voronoi dia- 
gram, seen from the negative y 
axis. 



'i-Z(W) o^ 

g(p) =\ 10 
1, 
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where \p\ is the parametric distance from p to the origin. The radial function / is non-negative, 
and it and its first two derivatives are bounded, e.g., 



f(r),\f'(r)\,\f"(r)\<0. 



(21) 



We also demand that there exists a positive 7 < j3 such that f(r) > 7 when r < ro, and that 
f(r) = if r > 2r . The parameter r , defines the radius of the ball bounding the perturbed 
region. Now we have d(w,p) < d(u,v) when £ < ro7. 

Since 7 may be arbitrarily small compared to /3, standard arguments supply a function / 
meeting these conditions. For example, the C°° construction described by Munkres |Mun68, p. 
6] may be multiplied by a scalar sufficiently small to meet our needs. 

The vertical y = cross-section of the perturbed Voronoi diagram will look something like 
Figure [3] Vj^{p) and Vj^(w) now meet in the nz-plane, and V^(u) and Vj^iv) do not. However, 
since geodesies which do not intersect the ball B^3 (0, 2r ) will remain straight lines in the pa- 
rameter space, the Voronoi diagram is unchanged outside of a neighbourhood of the origin. Thus 
looking from above at the slice of the Voronoi diagram in the xy-plane, we will see something 
like Figure ^[a) Figure |^[b) | shows the yz-plane. 

Two Voronoi vertices have been introduced, the red and 
blue points in Figure [4] These are the centres of dis- 
tinct empty geodesic circumballs for {p, u,v,w}. Since they 
cannot lie in the region unaffected by the perturbation, a 
quick calculation shows that the parametric distance of these 
Voronoi vertices from the origin is bounded by 4ro, when 
r o < 1 1 and it follows from another small calculation that 
the parametric distance from these Voronoi vertices to any of 





u 














• w 




• 






V 





The dis- 



Figure 3: The y = slice of the 
perturbed Voronoi diagram. 



the four sample points is bounded by a(l- 
tances between these Voronoi vertices and the sample points 
in the new metric will also be subjected to the same bound, 
since no distances increase. Also, The sparsity condition will 
not be affected by the perturbation. Thus, since we can make 
r as small as we please, and £ is chosen such that £ < r j, 
it follows that the radius of these balls may be made arbi- 



trarily close to a 



we can make 



feo?7(-M)- We will argue next that 



r](M)-r)(M) 



as small as desired by reducing 



the size of (3 in Equation (21 1. Then other sample points may be placed on the manifold so that 



the density criteria are met, and no degenerate configuration (violation of Definition A. 1 1 need 
be introduced. 

This means that the Delaunay complex, defined as the nerve of the Voronoi diagram, will not 
be a triangulation of the manifold M. . As observed by Boissonnat et al. |BGO09j , the triangle 
faces {p, w,u} and {p, w,v} will be adjacent to only a single tetrahedron, namely {p, u, v,w}. 
Thus De\^(V) is not a manifold complex as defined in Section [2] This is clearly a problem if 
the original manifold has no boundary. 

Although it is in some sense close to being degenerate, we emphasise that this configuration 
represents a problem that cannot be escaped by an arbitrarily small perturbation of the sample 
points. An argument based on the triangle inequality shows that in order to effect a change in 
the topology of the Voronoi diagram, a displacement of the points by a distance of f2(rc)7 ~ ^ s 
required. 

More specifically, we observe that the configuration {p, u, v, w} may be placed in an otherwise 
well behaved point set V such that within a small ball centred at the origin in our coordinate 
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(a) xy-plane from above 




Figure 4: Looking at cross-sections; the positive y-direction is to the right. The four points, 
p, u, v,w, admit two small circumballs with distinct centres (the red and blue points). 



chart, all points will have {p, u,v,w} as the four closest points in V, and this would remain the 
case even if the point positions were perturbed a small amount. We may further assume that the 
other Delaunay simplices are well shaped, so that stability results |BDG12| can be used to argue 
that they cannot be destroyed with an arbitrarily small perturbation. Then we argue that in 
order to obtain a triangulation by a perturbation V — > V' , we must ensure that the Voronoi cell 
"Km ({p'i w '}) mus t vanish: the edge {p' , w'} will never be incident to any tetrahedron other than 
{p' , u' , v' ,w'}. Then an argument based on the triangle inequality shows that for a p-perturbation 
with p < r ° , there will be a point in ({p', w'}) within a distance of 2p of the origin. 



A. 1.2 The sizing function under perturbation 

We need to establish that the metric manipulation that we performed in order to construct the 
counter-example, does not have a dramatic effect on the sizing function rj(M). This follows from 
the fact that we have bounded g — g together with its first and second deriviatives. 

Since the sectional curvature may be described as a continuous function of g and it's first and 
second derivatives |dC92, pp. 56 & 93], the effect of our perturbation on the sectional curvatures 



can be made arbitrarily small by reducing (3 in Equation (21 1 



Since we started with a flat metric anyway, the bound K.q can be made arbitrarily small, and 



so the second term in Equation (20) will not be the smallest. We need to bound the change in 
the injectivity radius as well. 

This follows from results in the literature |Ehr74, Sak83 , which state that for a compact man- 
ifold, inj(jM) depends continuously on the metric and its first and second derivatives. Specifically, 
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Lemma A. 4 (Ehrlich) Let CtJt be the space of C 3 Riemannian metric structures g on a compact 
manifold M.. and endow 9Jt with the C 2 topology. The function g n- inj g (.A4) is continuous in 
this topology. 



This means that for any desired bound on n{M.) — r)(JA) 
bound. 

The construction of the counter-example is complete. 



there will be a j3 that will satisfy the 



A. 2 Discussion 

We have shown that for constructing a Delaunay triangulation for an arbitrary Riemannian 
manifold, a sampling density requirement is not sufficient in general. The solution we propose 
in the body of this paper, is to constrain the kind of sample sets that we consider. Another 
approach would be to constrain the kind of metrics that are assumed. However, even with a 
purely Euclidean metric, allowing configurations to be arbitrarily close to degeneracy means that 
arbitrarily poorly shaped simplices are to be expected. When the metric is no longer Euclidean, 
the "shape" of a simplex no longer has an obvious meaning, but the problems associated with 
point configurations near degeneracy will certainly be present. 

Our analysis relied on the ability to make the support of the perturbation small. This is 
unlikely to be a necessary feature of the construction, but it facilitates our simplistic analysis. 

Clarkson [Cla06j remarked that an implication of Leibon and Letscher's claim [LLOO is that 
for four points close enough together, there is a unique circumsphere with small radius. Our 
counter-example shows that circumcentres need not be unique under these conditions. In fact 
the existence of unique circumcentres does not follow from the triangulation result: In our work 
we do not claim that the m-simplices have a unique circumcentre in the intrinsic metric. However, 
the argument sketched out by Leibon and Letscher claimed that the intrinsic Voronoi diagram 
is a cell complex (i.e., it satisfies the closed ball property [ES97 ), and this does imply unique 
circumcentres for the top dimensional simplices. 

It is worth emphasising that the problems discussed here only arise when the dimension 
is greater than 2. The same sampling criteria for two dimensional manifolds has been fully 
validated Lei99, DZM08], however these works both assume genericity in the sample set, without 
demonstrating that it is a reasonable assumption. 



B Background results for manifolds 

The tangent space at p £ M. is denoted T p M, and we identify it with an m-flat in the ambient 
space. The normal space, N p A4, is the orthogonal complement of T p M. in T p M. N , and we likewise 
treat it as the affine subspace of dimension m — k orthogonal to T p A4 C 1* . 

A ball B = B r n (c, r) is a medial ball at p if B n A4 = 0, it is tangent to A4 at p, and it 
is maximal in the sense that any ball which contains B either coincides with B or intersects 
A4. The local reach at p is the infimum of the radii of the medial balls at p, and the reach 
of A4, denoted rch(A4), is the infimum of the local reach over all points of M.. In order to 
approximate the geometry and topology with a simplical complex, manifolds with small reach 
require a higher sampling density than those with a larger reach. As is typical, an upper bound 
on our sampling radius will be proportional to rch(Af). Since M. C is a smooth, compact 
embedded submanifold, it has positive reach. 

An estimate of how the tangent space locally deviates from the manifold is given by an 
observation of Federer |Fed59| Theorem 4.8(7)] (see also Giesen and Wagner |GW04l Lemma 6]): 
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Lemma B.l (Distance to tangent space) If x,y € M. C R N and d K jv(x, y) < r < rch(A^), 
then ci R ,v (y, T X M) < 2l - c h(M) ' anc ^ thus sma — 2roiZM) ' wnere a ^ s the angle between [x,y] and 

A complementary result bounds the distance to the manifold from a point on a tangent space 
|BG10I Lemma 4.3]: 

Lemma B.2 (Distance to manifold) Suppose v € T X A4 with — a;|| = r 

< ££!ltMl. Let 

V = ^xiv) G AL where ^ x is the inverse projection ([6]). Then. d R N(v,y) < ^gfy^n • 

The previous two lemmas lead to a convenient bound on the angle between nearby tangent 
spaces. We prove here a variation on previous results [NSW08, Prop. 6.2] |BG11[ Lemma 5.5]: 

Lemma B.3 (Tangent space variation) Let x,y € M. be such that d R N(x,y) = r < rcll (^ /| ) 
and let a be the angle between T X A4 and T y A4. Then, sin a < 



6r 



4 



Proof Let u G T y M C M w with || 

v ~ y\\ = r - We will bound the angle between v — y and 

T X M. We have 



sina< - — - — - (d R N(y,T x M) + d R N(v,T x M)) 

\\ v -y\\ 

< j, M (d.RN(y,T x M) + d R N(v,v) + d R N(v,T x M)) , 

\\v-V\\ 



(22) 



where v € M. i s the closest point to v in AL 

we have d R N (y, T X M) < 2 rch(M) > anc ^ ^y Lemma B.2 we get d R N(v,v) < 



By Lemma 



B.l 



2r- 



ch(M) 



. For the third term in Equation (22 1, we find 



d R N(x,ii) < d R N(x,y) + \\v - y\\ +d R N(v,v) 

2r 2 5r 

< 2r -I < 

" rch(Af) - 2 



and so we may apply Lemma 



B.l 



< rch(Al), 
to obtain da 



Putting these observations back into Equation (22 1 we find 

1 / r 2 2r 2 25r 2 \ _ 45r 6r 

Sma - \\v-y\\ \2xcMM) + rch(M) + 8rch(A4) J ~ 8rch(Al) < rch(M) ' 

□ 

The following observation is a direct consequence of results established by Niyogi et al. 
|NSW08| Lemma 5.4]: 

Lemma B.4 Let W = B R N\ M (p,r). for some p g AJ and r < rch(A4)/2. When restricted to 
W, the orthogonal projection ir p \w '■ W — > T p A4 is a diffeomorphism onto its image. 

Proof Let / = 7t p |vk- Niyogi et al. showed [NSW08J Lemma 5.4] that the Jacobian of / is 
nonsingular on W, so that W is a covering space for U = f(W) C T p M.. The Morse-theory 
argument of Boissonnat and Chazals [BC01, Proposition 12 ] can be applied to demonstrate that 
W is a topological ball. It follows that U is connected, since any path in W projects to a path 
in U. Thus W must be a single-sheeted cover of U, since / _1 (0) = {p}. Indeed, if q € W with 



q ^ p and f{q) = 0, then [p, q] would be perpendicular to T p A4, contradicting Lemma B.l Thus 



/ : W — > U is a diffeomorphism. □ 



Inria 



Constructing intrinsic DTs 



45 



Niyogi et al [NSW08, Prop 6.3] demonstrate a bound on the geodesic distance between nearby 
points, with respect to the ambient distance. We will use a modified statement of this result: 



Lemma B.5 (Geodesic distance bound) Let x, y € M be such that d R j\r(x, y) < rch ^ A/1 ) 
Then 



d M (x,y)< d R N {x,y)(l+ ^/Xif 
Proof The announced result states 

under the same hypothesis on x and y. Rearranging, we have 

A , 2d R N(x,y) d MN (x,y) ( 2d m (x,y)\ 

dM{x > v) - : , l - i - dRN{x > y) { 1 + ^mmT ' 

1 + V vch(M) *MM) 

where the second inequality is obtained by squaring away the radical. □ 

C Forbidden volume calculation 

In this appendix we demonstrate: 

Lemma |5.8| (Volume of forbidden region) Let a be a /c-simplex with vertices on M and 
k < m. If 

L r o ^ bTT' 

2 - 1 ^ min { 4 gri(^) ' and 



3. <5 2 <mm{r™ +1 ,i} 



4^rch(X)' 8/3 

m+1 l-l 
) 4-1' 



then 

voi(F( ( x,t))<^r i?( ( Tr, 

where -D depends on m and /?. 



We will use the following lemmas in the proof of Lemma 5 



Lemma C.l (Triangle altitude bound) For any non-degenerate triangle a — [p,q,r], we 
have 

D{ P ,a)= ^- qUp - r K 
' 2R(a) 

Proof Let a — Zprq and observe that 

IIP-«Z1I 

Since D(p, a) = \\p — r\\ sin a, the result follows. □ 
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Lemma C.2 Let a — [po - • • Pfc] C K be a fc-simplex with 1 < k < m < N. Suppose 



Pk+i € R w is such that <Ji = Pk+i * o admits an elementary weight function uj ai : d\ — > [0, co), 



and the following conditions are satisfied: 

(1) L(a,) > |, 

(2) #Oi,w CT1 ) <0t, 

(3) o~i is a Fo-flake, and 

(4) 5 2 < min{r-+ 1 , i}. 
Then 

a R N ( Pk+1 ;dS')<BT R(a) 
where S' = B UN (C(a),R((r)) n aff(cr) and 

B d = 4 + 96/3(1 + 2 7 3 2 /3 2 ). 

Proof Let uj a — u) ai \&. Note that w CT : <r — > [0, oo) is an elementary weight function, and 
C(a,u! a ) is the o rtho gonal projection of C(cti,w CTi ) onto aff(cr). 

(2) and the fact that <5q < |, we have 



From Lemma 



4.1 



a(<ti) < YzrF R (ri> u "J < g-R^i.w^i). (23) 



and 



/ ' \ > ^ / o; V from Lemma O (2) 

i?(cri,o; (Tl ) 2i?(o-i,w ffl ) 1 1 

- sJ^! S as 6 o ^ I and ^ A ( CT ) 

8it(<Ti,iJ (Tl ) 4 

> i (24) 



24/3 

Fherefore, from Lemma 12.71 we have 



^+^<fl + l_W A( CT1 ) 2 



A(o-) V k J L(a 1 )A(a) 

< 2r x from fc > 1 and £(<n) < A(cr) 

L(cri)^ 

128r R(ai,Lo a ,) 2 „ „ ,— , 

< x from Eq. © 

< 2 7 3 2 /3 2 xr from hyp. [(1)1 fc [(2)1 (25) 

Let p be the point closest topfc + i in 8B^n(C; R) where C = C(ai,oj ai ) and R = R(ai,oj ai ). 
We have 

b-Pik+ill = VR 2 + ^ (Pk+i) 2 -R<Uai (Pfe+i) < <5oL(fTi) (26) 

Let q be the point closest to p on 9B r n (C; i?) n aff(cr), p' be the projection of p onto aff(cr), 
and and let r denotes the intersection of the line aS([qC(a,cj er )]) with 8B r n (C; R). Note that 
C(<ri, tJ ai ), C(a,uj a ), Pk+i, P, p' ' , q and r lie on the same 2-dimensional affine space. 
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Using the fact that \\p — Pk+i\\ < SqL{<j\), we get 

lb - Pi < £>(Pfc+ii o-i) + Wa x ) (27) 

We will now consider the triangle er 2 = [PQ r ]- Note that C(oi , uj ai ) , R(<7i, uj C t 1 ) are the 
circumcenter and radius of o~ 2 respectively. Also, C(a,u) a ) is the midpoint of the line segment 
[qr] with 2i?(er, ui a ) = \\q — r\\ and D(p,a 2 ) = \\p — p'\\- From the definition of q, we have 
\\p — r \\ > ||p — <?|| ■ Using the fact \\q — r\\ = 2i?(er, to a ), we have 



lb-r|| > L ^ = R{a,w <J ). 



This implies from Lemma |C.1| 
2R{a 2 )D(p, 0-2) 



\\p-q\\ 



\\p-r\\ 

< 2J?(gl D ^ l)g(P ' g2) as R(a 2 ) < Riat^i) and ||p - r|| > flfow,) 

H(a, uj a ) 

< 48(3D(p, a 2 ) = 48/% as Eq. p4l) (28) 



From Eq. (|26j), (|27j) and ([28 1 

~ q\\ < \\pk+i - p\\ + \\p - q\\ 

< ioL(ai) + 48/3(£>(p fc+1)( T 1 ) + * i(<n)) 

(29) 

Using the fact that a is Tp-thick (since Ui is a r -flake), and the bound 8qL((Ti) 2 on the 
differences of the squared distances between C(a, uj a ) and the vertices of er, we obtain a bound 
[BDG12I Lemma 4.1] on the distance from C(a, cj ct ) to C(a): 

\\C{a)-C{a^)\\< 5lL{ai)2 



T(cr) 

< *o#W since a is T^-thick, T(<r) > rg 

< v ; as r < 1 

1 

dcf 



»fc (30) 

Since fc > 1, there exists Pi <E a such that 

Pi G B H «-(C( C 7), J R(c7))nBH«(^(o-,w CT ), J R(c7,w CT ))naff((T). 

Also, ||C(cr) = R(a) and ||C(ct, - = R(a,LJ a ). 

Using the facts that R(a) = ||C(cr) — p»|| and R(a,uj a ) — \\C(a,u) a ) — Pi\\, and the Triangle 
inequality, we get 

R(tr)-\\C{*)-C(<T,w a )\\ < WCfawJ-piW < R(a) + \\C(a)-C{a,iJ a )\\ 

R(a)-r) 2 < R{a,uj a ) < R{a)+r ]2 (31) 
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The last equation follows from Eq. (|30|. 

Let S' and S denote _B R w (C(cr), -R(<j))naff (<r) and B^n (C(a, u> a ),R(a, u> a )) respectively. From 
Eq. (30) and (31 1, we have d R N (dS', dS) < r\\ +2rj2- This implies that there exists q' £ dS' such 
that 

y - q\\ < 2772- (32) 



Therefore from Eq. ( 29 ) and (32 I. we get 

\\Pk+i - < |bfe+i - q\\ 



Using the facts that 5$ < < (from hyp. [(4)] of the lemma and T < 1), L(oi) < 

L(a) < A(cr) < 2R(a) and D(p ^j <Tl) < 2 7 3 2 /3 2 T (from Eq. ([25)), and Eq. @ and Q, we 
get 



dR«0/c+i;<9S") < Ibk+i-g'l 

< r?i + 2r/ 2 



where 



< faLfa) + 48/3 (X>(pfc+i, ^l) + SoLfa)) + 

< BT Q R(a) 

B = 4 + 96/3(1 + 2 7 3 2 /3 2 ). 



2<5 2 fl(a) 



□ 



We will use the following lemma from [BG10 to bound the volume of F(a) 



Lemma C.3 Let p be a point on M.. There exists £ that depends on rch(TW) and m, and A 
that depends only on m such that, for all r = t < £, we have 

. ~ voI(B r n(p, r) f]M) 
< 1 - At < — ^ K K1 \' '- < 1 + At 

— V r k ~ 

where V m is the volume of the m-dimensional unit Euclidean ball. 

Proof of Lemma \5.8\ For the rest of the proof we define 

t 

t 



Tch(M) 

Consider the following elementary weight function: co a 



0J a 



Using the facts that 



R(o-,uj<j) < R{(Ti , uj ai ) , R(a,uj a ) < f3trch(Ai), and Lemma 4.1 (3) 

S 2 x _1 



R{a) < R(<j,w a ) 1- 



T(a) 



since T(o-) > r£ > 



< 2/3trch(7W) 



(33) 



Let p be a vertex of a. Let c be the point closest to C(a) on T V M. and c* be the point closest 
to c on M. (see Fig. [5]). 
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T P M, 




Figure 5: Proof of Lemma 5.8 



From Lemma |B.1| we have for all q e a 



A i T KA\S HP"*?!! 2 . A(CT) 2 dof 



From Lemma 2.1 and the facts that T(er) > T™ and A(cr) < 2R(<j) < 4/3i, we have 



s ri^aw s TOOT) s f 



Therefore 



Hc-CHH < sinZ(T p .M,ah») x i?(a) < i?(a), 



and from Lemma 14.21 



!<•-«■ II < 2||C , ?t < A0tR{&). 



(34) 



(35) 



rch(M) 

Let x € F(a,t) and a;* be the point closest to x on d B r n (C (a) , R(a)) D aff(cr). Then from 
Lemma [C.2[ we have 

||ic-a:*|| < BT Q R(a) (36) 
Using the fact that ||C(cr) — x*\\ — R(cr), we get 

||c* - x\\ < \\c* - c|| + ||c - C(a)\\ + \\C{a) ~x*\\ + \\x* - x\\ 



< R(a) ^1 + BT a + i/3t 

<*(-) ( 1 + ^+rf ) 
<i?(<7)(i + (s + i)r ) 

Similarly we can show that 



1 

J- n 



l 



from Eq. (34), (35), (36 1 



since r < 1 
from hyp. [2] of the lemma. 



\c*-x\\ < R(a)(l - (B + 1)T ) 
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Therefore 

F(a, t) C (B R , (c*, (1 + C)-R(o-)) \ (c*, (1 - W)) n A-f 

where £ = (B + 1)T Q . 

Observe that Lemma |C.3| can be applied since 

R(a)(l + C) < 2R{a) since C < 1 from hyp. [T] 

< Apt from Eq. ([33) 

< £ because i < e. 

Therefore 

vol(^M)) < vol(B R iv(c*, -R(g)(l + Q) n M \ -Brw (c*, fl(<x)(l - Q) n M) 

< (1 + A(l + C)t)R{<j) m {l + C) m - (1 - A(l - ()i)R(a) m (l + cr 

< R{<jy n {(i + cr - (1 - cd + AtR{<j) m {{\ + o m + (1 - en 

< 2 m ( R{a) m + A{2 rn+1 + l)tR{a) m (37) 

The last inequality follows from the fact that (1 + x) m - (1 - x) m < 2 m x for all x € [0, 1]. 
From hyp. [2] and the fact that T < 1, we have 

pm+l 

t<l<^ r <T Q . (38) 



The lemma now follows from Eq. ( 37 1 and ( 38 ) . □ 
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